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Preface 



The aim of this book is to deal with all of the elementary methods for obtaining 
explicit solutions of ordinary differential equations, and then to introduce the ideas 
of qualitative analysis using phase plane techniques. Simple difference equations 
are also included, since their methods of solution are similar to those for linear 
differential equations. As well as being, I hope, an internally consistent choice of 
material, this selection of topics also has the advantage of preparing a student for 
a basic course on dynamical systems. 

The book arose from my unsuccessful efforts to find a suitable text to recom- 
mend when I taught the first year Warwick differential equations course. Although 
there are a number of well-established and successful textbooks that treat this sub- 
ject (these are discussed, along with other possibilities for further reading, in the 
final chapter), they seem either to include a large amount of additional material, or 
to concentrate only on the more advanced topics. I therefore produced a detailed set 
of lecture notes, which, with the encouragement of Alan Harvey and David Tranah, 
and most significantly Kenneth Blake at Cambridge University Press, eventually 
became this book. My thanks here to all those students who made useful sugges- 
tions while this book was still at the lecture note stage. 

Part I contains an informal discussion of the issues of existence and uniqueness 
of solutions, and treats the standard classes of first order differential equations 
that can be solved explicitly, as well as covering exact equations and substitution 
methods. 

The first chapter of Part II shows that two linearly independent solutions are 
needed in order to solve the general homogeneous problem, and also contains a 
brief treatment of the Wronskian. The remainder of this section treats equations 
with constant coefficients, concentrating for the most part on the second order 
case, with higher order equations discussed briefly at the end. 

Second order equations with non-constant coefficients are treated in Part III, 
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xiv Preface 

which covers reduction of order, the method of variation of constants, and series 
solutions. 

Part IV turns aside from differential equations, motivating the study of dif- 
ference equations by discussing Euler’s method of numerical solution. Constant 
coefficient linear difference equations are covered, and then there are two chapters 
devoted to nonlinear difference equations. One of these goes beyond the coniines 
of an introductory course and discusses the dynamics of the logistic map in some 
detail. 

Part V treats coupled systems of two linear differential equations, starting with 
the substitution method that reduces the problem to a second order differential 
equation in one variable, the most reliable way to find explicit solutions. The re- 
mainder of this portion of the book deals with the matrix approach, showing how a 
calculation of the eigenvalues and eigenvectors of an appropriate matrix is enough 
to draw the phase portrait. This is done by changing to a coordinate system in 
which the equation is put into a standard form, providing an illustration of the 
Jordan canonical form of a matrix. 

Part VI uses the methods from Part V in order to draw the phase plane diagrams 
for a variety of nonlinear systems, with examples taken from mathematical ecology 
and simple one-dimensional particle systems, including the pendulum. The book 
ends with a brief discussion of Dulac’s criterion and the Poincare-Bendixson The- 
orem, a chapter that investigates the complicated dynamics of the Lorenz Equa- 
tions, and suggestions for further reading. 

In addition to those already mentioned above I would like to thank various peo- 
ple who have contributed to this book. I first learned much of the material here 
from Tristram Jones-Parry at Westminster School, to whom much belated thanks 
for all his fine teaching many years ago. I also owe a debt of gratitude to all those 
who taught the course at Warwick before me, shaping its contents and therefore 
those of this book; in particular, I had useful guidance from the course notes of 
Alan Newell and Claude Baesens. I am most grateful to Andrew Stuart, who, in 
encouraging me to emphasise the links with linear algebra, made me fond of a 
subject that I still remembered with a shudder from my own undergraduate days. 
Thanks too to James Macdonald, whose ‘Swarm of flies’ program for his MMath 
project on the Lorenz equations was the inspiration behind Figure 37.8. 

Over the past two months I have been able to think of little except phase planes 
and drawing figures in Matlab: my wife, Tania Styles, has managed to endure my 
many variations on ‘come and see this picture of a washing machine’ with a smile. 
Heartfelt thanks to her for this, and, of course, for everything. 

Finally, I would particularly like to thank my Ph.D. student, Oliver Teame, and 
my father, John, both of whom read this book extremely carefully and made a num- 
ber of very helpful comments. For whatever imperfections remain, my apologies 
to them and to my readers. 




Introduction 



Differential equations date back to the mid-seventeenth century, when calculus 
was discovered independently by Newton (c. 1665) and Leibniz (c. 1684). Mod- 
em mathematical physics essentially started with Newton’s Principia (published in 
1687) in which he not only developed the calculus but also presented his three fun- 
damental laws of motion that have made the mathematical modelling of physical 
phenomena possible. 1 

Historically, advances in the theory of differential equations have come from the 
insights gained when trying to treat specific physical models. Despite this some- 
what piecemeal development, the subject has become a well-defined and coherent 
area of mathematics. This book adopts a theoretical point of view, developing the 
theory to the point at which it can no longer be described as ‘basic differential 
equations’ and is about to become entangled with more advanced topics from the 
theory of dynamical systems. Of course, applications are used throughout to serve 
as motivation and illustration, but the emphasis is on a clean presentation of the 
mathematics. 

You may find that some of the problems covered in the first few chapters are 
already familiar. The methods of solving these problems are well established, and 
you may be well practised at applying them. However, we will take care here to 
show why these methods work; giving proper justification of the methods can take 
some time, but as mathematicians we should not be satisfied merely with a set of 
‘recipes’. Nevertheless, knowing something about the details should not stop you 
from applying the methods you know already; rather you should be able to use 
them with more confidence. 

Some of the chapters, and some sections within other chapters, are marked 
with an asterisk (*). These parts of the book contain either material that is more 
advanced, or material that expands on points raised elsewhere; while they could be 
omitted in the interests of brevity, they are intended to give some indication of the 
richness of the subject beyond the confines of an introductory course. 

1 Various modem editions of this work are available, translated from its original Latin. 
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Introduction 



There are three appendices, covering background material that is necessary at 
various points in the book. While some of this is elementary and may already 
be familiar (Appendix A recalls some notation and various facts about real and 
complex numbers that will be used throughout the book) some is a little more ad- 
vanced. Problems with timetabling often mean that certain undergraduate courses 
have to rely on material that is yet to be taught in others, hence there are appen- 
dices on matrices, eigenvalues and eigenvectors (Appendix B) and on derivatives, 
partial derivatives and Taylor series (Appendix C). The calculation of eigenvalues 
and eigenvectors is treated in detail in the main part of the book. 

The use of mathematical computer packages is now a standard part of the under- 
graduate curriculum, and an important tool in the armoury of practising mathemati- 
cians, scientists and engineers. Although the emphasis in the text is on pencil and 
paper analysis, and the book in no way relies on the availability of such software, 
some topics, particularly the treatment of coupled nonlinear equations using phase 
plane ideas in Chapters 28-37, can benefit greatly from the graphical possibilities 
modem computers provide. Almost all of the figures in this book have been gener- 
ated using Matlab, and very occasionally particular Matlab commands are men- 
tioned in the text. Nevertheless, it should be possible to carry out the numerical ex- 
ercises suggested here using any of the major commercially available mathematical 
packages; and with a little more ingenuity using any programming language with 
graphical capabilities. The Matlab files used to produce some of the figures, and 
mentioned in certain of the exercises, are available for download from the web at 
www. Cambridge . org/ 0521533910 . 

There is no better way to leam this material than by working through a selection 
of examples. One set of examples is included in what is, I hope, a natural way in 
the text, with the end of each worked solution marked with a box (□). Another set 
of examples is given in the exercises that end each chapter, and these should be 
considered an integral part of the book. The majority consist of sample problems 
that can be treated with the methods of the chapter - in order to give teachers a 
reasonable choice of problems, there are intentionally more of these than you could 
reasonably be expected to do. Others, labelled with a ‘T’, are more theoretical 
and designed to give an indication of some of the mathematical issues raised, but 
not treated in detail, in the text. Finally, those exercises labelled with a ‘C’ are 
intended to encourage the use of the computer to perform routine calculations and 
investigate equations and their solutions graphically. Those involved in teaching 
courses based on this book may obtain copies of solutions to these exercises by 
applying to the publisher by email (solutionsOcambridge . org). 

I would welcome any comments or suggestions, either by post to the Mathe- 
matics Institute, University of Warwick, Coventry, CV4 7AL, U.K. or by email 
to j crSmaths . Warwick . ac . uk; any errata that arise will be posted on my own 
website www. maths . Warwick . ac . uk/~ j cr/IntroODEs . html. 
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First order differential equations 
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Radioactive decay and carbon dating 



Before we start our formal treatment of the subject we will look at a very simple 
example that nonetheless exhibits the power of differential equations as models of 
reality. One point to bear in mind in this chapter is the distinction to be made be- 
tween finding the solution of a differential equation, and interpreting this solution. 



1.1 Radioactive decay 

Let N(t) denote the number of radioactive atoms in some sample of material at 
time t. Then with k > 0 the equation 

dV 

— = ~kN (1.1) 

at 

is a very good model for the way that the number of radioactive atoms decays (see 
Exercise 1.1). 

Although we will see later how to solve this equation, for now we will assume 
that when there are N s isotopes at time s, the solution is 

N(t) = N s e- k(t - s) . (1.2) 

You can check that we really do have the solution: when t = s the formula in (1.2) 
gives N ( 5 ) = N s , while we have 

—N(t) = -kN,e~ k(t ~ s) = —kN(t), 
at 

and so the differential equation (1.1) is satisfied. 

It follows from (1.2) that the number of radioactive isotopes decays exponen- 
tially to zero. Graphs of the solution for various values of /V (0), showing this 
decay, are plotted in Figure 1.1. 

The half-life of a particular radioactive isotope is the time it takes for half of the 
radioactive isotopes to decay, and this is related to the constant k that appeal's in 
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1 Radioactive decay and carbon dating 




Fig. 1.1. Graph showing the number N(t ) of radioactive atoms falling off as a 
function of time, for a number of different values of /Vo; the constant k is that for 
radioactive carbon 14. The half-life, approximately 5700 years, is marked by a 
dashed vertical line. 

the equation. To find this relationship, suppose that there are No radioactive atoms 
at time t = 0. Then the solution of (1. 1) is 

N(t) = N 0 e~ kt . 

Half of the atoms will have decayed by time thaif when N (/half) = o A/), i.e. 

N 0 e~ ktialf = j A/) =4> e-^ h aif = 1. 

Taking the (natural) logarithm of both sides gives 

-fahalf = -In 2, 

and so the half-life is given by f half = (In 2)/ k. Note that this time does not depend 
on the initial number of radioactive atoms. 



1.2 Radiocarbon dating 

The solution (1.2) forms the basis of the technology of radiocarbon dating. The 
essence of the method is as follows. Living matter is constantly taking up carbon 
from the air. The result is that within such material the ratio of the number of iso- 
topes of radioactive carbon 14 ( 14 C) to the number of isotopes of stable carbon 12 
( 12 C) is essentially constant. Once the specimen is dead (for example, a tree is cut 
down for its wood, or cotton is harvested for weaving), the radioactive 14 C atoms 
begin to decay according to the model (1.1). Since the half-life of carbon 14 is 
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approximately 5700 years, we need to take the constant k in (1. 1) to be 

In 2 . 

k = as 1.216 x 10“ 4 . 

5700 

By examining the ratio of the number of isotopes of carbon 12 to carbon 14 
in a sample of the material that we want to date, it is possible to work out the 
proportion remaining of the 14 C atoms that were initially present. Suppose that the 
sample stopped taking up carbon from the air when t = s, and that the number 
of 14 C atoms present then was N s . If we know that the sample now (at time to) 
contains only a fraction p of the initial level of 14 C, then N (to) = pN s . 

Using our explicit solution N(t) = N s e~ k(t ~ s \ we should have 

pN s = IV (t 0 ) = N s e~ k{t °- S \ 

Cancelling the factor of N s in the two outside terms yields the equation 

p = Q- k ^- s \ 

Taking logarithms of both sides we have 

In p = —k(to — s), 



and so the year s from which the sample dates is given by 

In p 

s = to H — - — • 
k 



(1.3) 



In 1988, the Shroud of Turin (see Figure 1.2) was dated by three independent 
groups of scientists from Arizona, Oxford and Zurich. Fibres from the shroud were 




Fig. 1.2. The Shroud of Turin: carbon dated to the fourteenth century. Photograph 
© 1978 Barrie M. Schwortz (his website at www . shroud . com is well worth a 
visit). 
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1 Radioactive decay and carbon dating 



found to contain about 92% of the level in living matter. 1 Using the expression in 
(1.3) shows that the Shroud therefore dates from 



5 = 1988 + 



In 0.92 

0.000 121 6 






1302, 



putting its origin squarely in the Middle Ages. 



Exercises 

1 . 1 Radioactive isotopes decay at random, with a fixed probability of decay per unit time. 
Over a time interval At, suppose that the probability of any one isotope decaying is 
k At. If there are N isotopes, how many will decay on average over a time interval At? 
Deduce that 

N(t + At) - N(t) « —NkAt, 

and hence that dN/dt = —kN is an appropriate model for radioactive decay. 

1 .2 Plutonium 239, virtually non-existent in nature, is one of the radioactive materials used 
in the production of nuclear weapons, and is a by-product of the generation of power 
in a nuclear reactor. Its half-life is approximately 24 000 years. What is the value of k 
that should be used in (1.1) for this isotope? 

1.3 In 1947 a large collection of papyrus scrolls, including the oldest known manuscript 
version of portions of the Old Testament, was found in a cave near the Dead Sea; they 
have come to be known as the ‘Dead Sea Scrolls’. The scroll containing the book of 
Isaiah was dated in 1994 using the radiocarbon technique; 2 it was found to contain 
between 75% and 77% of the initial level of carbon 14. Between which dates was the 
scroll written? 

1.4 A large round table hangs on the wall of the castle in Winchester. Many would like 
to believe that this is the Round Table of King Arthur, who (so legend would have it) 
was at the height of his powers in about AD 500. If the table dates from this time, 
what proportion of the original carbon 14 would remain? In 1976 the table was dated 
using the radiocarbon technique, and 91.6% of the original quantity of carbon 14 was 
found. 3 From when does the table date? 

1.5 Radiocarbon dating is an extremely delicate process. Suppose that the percentage of 
carbon 14 remaining is known to lie in the range 0.99 p to 1 .01 p. What is the range of 
possible dates for the sample? 

1 P. E. Damon et at . 'Radiocarbon dating of the Shroud of Turin' , Nature 337 ( 1 989), 6 1 1-6 1 5 . 

2 A. J. Jull et at.. 'Radiocarbon dating of the scrolls and linen fragments from the Judean Desert’, Radiocarbon 
37 (1995), 11-19. 

3 M. Biddle, King Arthur’s Round Table (Boydell Press, 2001). 
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Integration variables 



Because of the intimate relationship between differentiation and integration 
(discussed in more detail in the next chapter) there will be many integrals in this 
book, and it is worth pausing now in order to make sure that we have an appropri- 
ately unambiguous notation. 

Although in theory mathematicians make careful distinctions between ‘the func- 
tion /’ and ‘ /(*)’, the value that / takes at a particular point x, this distinction is 
rarely maintained in day-to-day informal discussions. 

Usually this does not cause any trouble. However, consider the following prob- 
lem, posed in ‘everyday’ language: 

Find the area under the graph of f(x) between a and x. 

Although the meaning of this is clear, ‘find the shaded area in Figure 2.1’, there is 
some potential for confusion when we try to write this down mathematically, since 
there arc too many vs around. Converting the English into symbols gives 

fix) dv, (2.1) 

and it should be clear that this is not satisfactory, since the symbol v is used in two 
different ways: once as the upper limit of the range of integration ( f x ), and once 
as the variable that is being integrated over (dv). 

When we integrate a function between two limits, for example 1 

f fix) dx, 

J a 

the variable that we are integrating over is a ‘dummy’ variable. It is just there to 
tell us how to do the integration, and plays no role in the final answer, which will 

1 Observe that there is no need to change our notation for this particular definite integral, since no confusion can 
arise as to the role of x. 
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2 Integration variables 




Fig. 2.1. ‘Find the shaded area’. 



only depend on a and b. So 

f b fix) dx = f f(0) de = [ b /(H) dH. 

J a J a J a 

(We can change the name of the dummy variable with no effect on the integral.) 

The obvious solution, then, is to change the integration variable in (2. 1) to some- 
thing other than x. However, changing the variable to something completely dif- 
ferent from v is likely to be confusing. The approach we will adopt will be to add 
a tilde ~ to the integration variable, so that instead of (2. 1) we will write 




fix) dx. 



( 2 . 2 ) 



All being well this should keep things ‘clean’ but should not be too jailing. 

We will also do something similar when evaluating integrals where v is an upper 
limit, i.e. 




fix) dx = 



Fix) 




when F' = f. 

Of course, very few people are this careful when they are doing calculations and 
the backs of mathematicians’ envelopes arc full of things like (2.1) rather than the 
pedantic (2.2). 





3 

Classification of differential equations 



Before we begin we need to introduce a simple classification of differential equa- 
tions which will let us increase the complexity of the problems we consider in a 
systematic way. 

3.1 Ordinary and partial differential equations 

The most significant distinction is between ordinary and partial differential equa- 
tions, and this depends on whether ordinary or partial derivatives occur. 

Partial derivatives cannot occur when there is only one independent variable. 
The independent variables arc usually the arguments of the function that we arc 
trying to find, e.g. x in fix), t in x(t), both v and y in G{x , y). The most common 
independent variables we will use are x and t, and we will adopt a special short- 
hand for derivatives with respect to these variables: we will use a dot for d/d t, so 
that 




and z = 



d ^ 

d t~ ' 



and a prime symbol for d/clt, so that 



y 



t 



dy 

dx 



and 



d 2 y 

dx 2 



Usually we will prefer to use time as the independent variable. 

In an ordinary differential equation (ODE) there is only one independent vari- 
able, for example the variable v in the equation 



dy 

dx 



= fix), 
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3 Classification of differential equations 



specifying the slope of the graph of the function y; the variable t in 

mx = f (t) 



which we could solve for the position x(f) = ( x(t ), y{t),z(t)) of a particle at 
time t moving under the action of a force f(t) (the equation is Newton’s second 
law of motion, F = ma); or x in 

h 2 d 2 xfi 

- — -fi r + V(x)fi = Efi 
2m dx - 

where i// (x ) = u(x) + i ft (x ) is complex (this is the Schrodinger equation from 
quantum mechanics). 

In a partial differential equation there is more than one independent variable 
and the derivatives are therefore partial derivatives, for example the heat in a rod 
at position x and time t,h{x,t), obeys the heat equation 

3/7 _ d 2 h 
3 1 dx 2 

A much more complicated example is given by the Navier-Stokes equations used 
to determine the velocity of a fluid 



u(xi, X2, X3, t) = (ll \(x\, X2, X3, t ), U 2 (xi, X2, X3, t ), ll 3 (x\ , X2, X3, t )) 
(think of xi = x, X2 — y, and X3 = z), which are: 1 



P 



3 iff 
3 1 



+ E 



Uj- 



\i = 1 



3 uj_ 

3 xi 



ix 



d~u j 
~dxf 



+ 



d-u , 



dx 2 



+ 



9 J*± 

dx 2 



(one for each component, j = 1,2,3) and 

du\ din dux 

— + — + — = 0. 

3xi 3x2 dx^ 



+ Y~ = fj- (3-D 

3x ; 



(3.2) 



In this book we will consider only ordinary differential equations. 



1 It is possible to write these two equations much more concisely using vector calculus notation. Imagine that V 

represents a vector of partial derivatives, V = (3/3 a'i , 3/3x2, 9/9x3), which can be manipulated like a normal 
vector. Then, for example. Equation (3.2) is just V ■ u = 0. Defining also A = V • V (the sum of all second 
derivatives) we can rewrite (3.1) as 



“ 3 u 

fit 



+ 



(u ■ V)u 



P 



- n Au + vp = f. 



(3.3) 




3.3 Linear and nonlinear 
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3.2 The order of a differential equation 

The order of a differential equation is the highest order derivative that occurs: the 
equation 



dy 

dx 



= fix) 



specifying the slope of a graph is first order, as is the following equation expressing 
energy conservation. 



\mx 2 + V (x) = E 



(\mx 2 is the kinetic energy while V (x) is the potential energy at a point x); 
Newton’s second law of motion 




F 



is second order; the equation 



y + \xfrxj/' = 0 



(which occurs in the theory of fluid boundary layers) is third order (recall that x//'" 
is shorthand for dS///dx 3 ). 

To be more formal, an nth order ordinary differential equation for a function 
y(t) is an equation of the form 



(d n _y d n _^y d y \ =0 

\dt" ’ dt' ,_1 ’ ’ d t ) 



(3.4) 



(Of course we want d n y/dt n to occur in 3F\ if Tiy, y, y, t) is y — t then the result- 
ing equation (y — t = 0) is not a differential equation at all.) If t does not occur 
explicitly in the equation, as in 



dy 

df 



= fiy), 



then the equation is said to be autonomous. 



3.3 Linear and nonlinear 

Another important concept in the classification of differential equations is lineality. 
Generally, linear problems are relatively ‘easy’ (which means that we can find an 
explicit solution) and nonlinear problems are ‘hard’ (which means that we cannot 
solve them explicitly except in very particular cases). 




14 3 Classification of differential equations 

An nth order ODE for y(t) is said to be linear if it can be written in the form 

d"y d' ,_1 y dv 

a n(t) — + a n -\{t)——r H \-ai(t)— +a 0 (t)y = /(?), (3.5) 

at" d t n ~ v at 

i.e. only multiples of y and its derivatives occur. Such a linear equation is called 
homogeneous if fit) = 0. and inhomogeneous if fit) f 0. 



3.4 Different types of solution 

When we try to solve a differential equation we may obtain various possible types 
of solution, depending on the equation. Ideally, perhaps, we would find a fully ex- 
plicit solution , in which the dependent variable is given explicitly as a combination 
of elementary functions of the independent variable, as in 

y(t) — 3 cos 5? + 8sint. (3.6) 

We can expect to be able to find such a fully explicit solution only for a very limited 
set of examples. 

A little more likely is a solution in which y is still given directly as a function 
of t, but as an expression involving an integral, for example 

y(t) = 1 + [ a~ s2 ds. (3.7) 

Jo 

Here y is still an explicit function of t, but the integral cannot be evaluated in terms 
of elementary functions. 

Sometimes, however, we will only be able to obtain an implicit form of the 
solution; this is when we obtain an equation that involves no derivatives and relates 
the dependent and independent variables. 2 For example, the equation 

In y + 4 In v — y — 2x + 4 — 0 (3.8) 

relates x and y, but cannot be solved explicitly for y as a function of x. 

All these types of solution will occur in what follows. 

There are many situations, however, in which it is not possible to obtain any use- 
ful expression for the solution. For some equations it is still possible to understand 



- We could also have an implicit solution containing integrals that cannot be evaluated in terms of elementary 
functions. For example, we will see that the equation dx/dt = f(x)g(t ) has solution 



[*L=f 

J fix) J 



= / g(t)dt, 



which in general gives such an implicit solution. 




3.4 Different types of solution 
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Fig. 3.1. A qualitative, graphical solution of the coupled system of equations 
(3.9). The axes are x (horizontally) and y (vertically), and it is safe to assume 
that this is the case for any unlabelled axes in the rest of the book. 



the qualitative behaviour of the solutions, i.e. to describe how the solutions be- 
have, even though we cannot specify them exactly. This is the approach we will 
take in Chapter 7, and throughout Chapters 32-37. Such a description is often best 
expressed graphically. For example. Figure 3. 1 shows the phase diagram (or phase 
portrait ) for the solutions of the equations 

x = x(4-2x-y) (39) 

v - y(9 — 3.v — 3y). 

The diagram is a plot of sample curves traced out by solutions (x(t), y(t)) labelled 
with arrows indicating the direction in which t increases. The crosses show points 
at which the solutions of this equation are constant. We can tell from this dia- 
gram that every solution eventually approaches the point (1,2) [i.e. x(t) 1 and 
y(t) — 2 as t — > +oo], even though we do not have any form of explicit solution 
for (3.9). 

For some equations all our analytical tools may fail, and in this case we can 
often use a computer to approximate the solution. A ‘numerical solution’ of a 
differential equation is usually only an approximation, and the initial result of such 
a calculation will not be an expression for x in terms of t, say, but a list of times, 
t , and corresponding approximate values for x(t). Using Matlab’s ODE solving 
routine, ode45, to solve the equation 




x (0) = 0 
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3 Classification of differential equations 



between times t — 0 and t — 5, yields such a list: 
>> xdot=inline ( ' t-x~2 ' , ' t' , 'x' ) ; 



» [t 


x] =ode45 (xdot , [ 0 5 ] , 0 ) ; 






» [t 


x] 








ans = 


0 


0 








0.1250 


0 . 0078 


2 . 6250 


1.4921 




0.2500 


0 . 0312 


2.7500 


1 . 5407 




0.3750 


0 . 0700 


2 . 8750 


1.5864 




0.5000 


0.1235 


3 . 0000 


1.6299 




0 . 6250 


0.1907 


3 . 1250 


1.6721 




0.7500 


0.2700 


3.2500 


1 . 7127 




0.8750 


0.3591 


3.3750 


1 . 7515 




1.0000 


0.4555 


3.5000 


1.7891 




1.1250 


0.5563 


3 . 6250 


1.8261 




1.2500 


0 . 6585 


3.7500 


1.8621 




1.3750 


0.7596 


3 . 8750 


1.8969 




1.5000 


0.8574 


4 . 0000 


1.9310 




1 . 6250 


0.9505 


4 . 1250 


1.9646 




1.7500 


1 . 0377 


4.2500 


1.9976 




1.8750 


1.1187 


4.3750 


2 . 0297 




2 . 0000 


1.1935 


4.5000 


2 . 0612 




2.1250 


1.2628 


4 . 6250 


2 . 0925 




2.2500 


1.3268 


4.7500 


2.1231 




2.3750 


1.3856 


4 . 8750 


2.1531 




2.5000 


1.4403 


5.0000 


2.1826 



We will discuss one simple method of numerical approximation in Chapter 2 1 . 



Exercises 



3.1 Classify the following equations as ordinary or partial, give their order, and state 
whether they are linear or nonlinear. In each case identify the dependent and inde- 
pendent variables. 

(i) Bessel’s equation (v is a parameter) 

x 2 y" + xy' + (.r 2 — v 2 )y — 0, 



(ii) Burger’s equation (v is a parameter) 



3 u d 2 u 3 u 

v — + u — 

3 t dx 2 dx 



= 0 , 




Exercises 
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(iii) van der Pol’s equation ( m , k, a and b are parameters) 

mix + kx — ax — bx , 

(iv) dy/dt = t — y 2 , 

(v) the wave equation (c is a parameter) 

d 2 y 2 d 2 y 

3 1 2 dx 2 ’ 

(vi) Newton’s law of cooling ( k is a parameter and Alt) is a specified function) 

d7 

— = -k(T - A(t)), 
at 

(vii) the logistic population model (k is a parameter) 

dp 

— = kp(l - p), 

(viii) Newton’s second law for a particle of mass m moving in a potential V (x ) , 

mx = —V'(x), 

(ix) the coupled equations in (3.9) 

x — x(4 — 2x — y) 
y = y(9 — 3x — 3y), 

and 

(x) 



where x is an n -component vector and A is an n x n matrix. 




4 

^Graphical representation of solutions using Matlab 



The list of numbers that formed the example of a numerical solution at the end 
of the previous chapter indicates how useful a graphical representation of solu- 
tions can be. In fact Matlab ’s default presentation of a numerical solution of a 
differential equation is as a graph: the commands 

>> xdot=inline ( ' t-x~2 ' , ' t ' , 'x'); 

>> ode45(xdot, [0 5], 0) 

produce the graph shown in Figure 4. 1 (only the axis labels have been added). 

Whichever kind of solution we manage to obtain for our equation, the graph- 
ical capabilities provided by modern computer packages enable us to visualise 
these solutions and so obtain a much better understanding of their behaviour. All 
the solutions in Section 3.4 benefit from a graphical presentation. In this section 
we briefly discuss the main Matlab commands that can be used to visualise and 
solve a variety of equations. 

Almost all of the figures in Parts I, II, and III of this book are the graphs of 
explicit solutions; these arc very easy to produce with Matlab. For example, to 
plot y(t) = 3 cos 5? + 8 sin t against t for 0 < t < 20, the three lines 

t=linspace ( 0 , 2 0 ) ; 
y=3*cos(5*t)+8*sin(t) ; 
plot (t,y) 

produce Figure 4.2. 

If the solution is given as an integral that cannot be evaluated explicitly, like 
(3.7), 



y(t) = l + 




ds. 
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* Graphical representation of solutions using MATLAB 
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Fig. 4.1. The solution of x = t — x 2 with x(0) = 0, as produced by the 
MATLAB ode45 command. The individual pairs (x, t) are represented by the 
circles, and are joined to produce an approximation to the solution x(t) of 
the original equations. 




Fig. 4.2. The graph of y(t) — 3 cos 5t + 8 sin t (y against t). 

then we can find the value of y at any given value of t by approximating the in- 

2 

tegral; this is something that computers are very good at. The integral of e -r 
between 0 and 2 (for example) can be evaluated by defining an ‘inline function’ 
f{t) — exp (— f 2 ) and then using the quad command: 

>> f =inline ( ' exp ( -t . ~2 ) ' , ' t ' ) 



f = Inline function: 



f ( t ) = exp ( - 1 . ~ 2 ) 






Exercises 
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Given an implicit formula like (3.8), 

In y + 41nx — y — 2x = —4, 

we can notice that x and y lie on a curve that makes 

F(x, y ) = In v + 41nx — y — 2x 

constant. The ‘contour plot’ of the level set F(pc, y ) = —4, 

» [x, y] =meshgrid ( . 01 : . 1 : 5 , . 01:. 1:5); 

>> z=log (y) +4*log (x) -y-2*x; 

>> contour (x, y , z , [ -4 -4]) 

is shown in Figure 4.4. 

Exercises 



4.1 (C) Plot the graphs of the following functions: 

(i) y(t) = sin5f sin50f forO < t < 3, 

(ii) x(t) = e _f (cos2 1 ± sin 2t) forO < t < 5, 

(iii) 

T(t) = f sin.vd.s- for 0 < t < 7, 

Jo 

(iv) x(t) — t In t for 0 < t < 5, 

(v) plot y against x, where 

x(t) — Be~' + Ate~ f and y(t) = Ae~ l , 



for A and B taking integer values between —3 and 3. 
4.2 (C) Draw contour plots of the following functions: 

0 ) 



(ii) 



F(x, y) = x 2 + y 2 for — 2 < x, y < 2; 



F(x, y) — xy for — 1 < x, y < 1, 



with contour lines where F — ±0.1, ±0.2, ±0.4, and ±0.8; 



(iii) 



E(x, y) = y“ — 2cosx for — 4 < x, y < 4; 



(iv) 



E(x, y) — x — |x 3 ± \y 2 (x A — 2x 2 ± 2) 

for —2 < x < 4 and —2 < y < 2, showing contour lines where E = 0, 0.5, 
0.8, 1, 2, 3 and 4; 



E(x, y) — y 2 ± x 3 — x 



for 



— 2 < x, y < 2. 
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‘Trivial’ differential equations 



In this chapter we consider the simplest possible kind of differential equation, 
one that can be solved directly by integration. Although the problem is relatively 
straightforward, it will serve to introduce several important ideas. 

You have probably already met and solved one simple kind of differential equa- 
tion: 

= fix). (5.1) 

ax 

Viewed as an equation to solve for y(x), this asks us to find the function whose 
graph has slope f (x) at the point v. So in order to solve this equation we ‘just’ 
have to find a function whose derivative is fix). 



5.1 The Fundamental Theorem of Calculus 

Any function F that satisfies F' = f is called an anti-derivative 1 of /. Clearly if 
F is an anti-derivative of / then so is Fix) + c for any constant c. 

This terminology allows us to distinguish between reversing the process of dif- 
ferentiation (finding an anti-derivative) and integration (finding the area under a 
curve). Put like this it becomes possible, perhaps, to appreciate how remarkable 
it is that these two concepts are so intimately related. This is formalised in the 
Fundamental Theorem of Calculus (FTC). 

Essentially this theorem says that differentiation reverses the action of integra- 
tion, and that if we know an anti-derivative of / we can calculate the area under 
the graph of / between any two points; it is easy to forget that the FTC is a major 
result because we use it so frequently in order to calculate integrals. 

In the statement of the theorem we use M to denote the set of all real numbers, 
and [n, b\ denotes the closed interval a < x < b. 

1 The more puzzling word ‘primitive’ is sometimes used instead of ‘anti-derivative’. 
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5.1 The Fundamental Theorem of Calculus 
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Theorem 5.1 Suppose that f : [a , b] — »• R is continuous, and for a < x < b 
define 



G(x)= f f(x)dx (5.2) 



(the integral G (x ) is the area under the graph of f between a and x, see 
Figure 5. 1 ). Then 

d G 

— (x) = f(x). 



and furthermore 



f 



f(x) dx = Fib ) — F(a ) 



for any anti-derivative F of f (i.e. for any F with F' = f). 



(5.3) 



We often write (5.3) in the more convenient shorthand 



f 



fix) dx = 



Fix) 



(5.4) 



Proof (Sketch) If we calculate dG /Ax using the formal definition of the derivative 
as a limit (see Appendix C), 



dG 

dx 



(x) = 



lim 

8x— m 



G(.r + bx ) — G(x) 
hx 



(5.5) 
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5 ‘Trivial’ differential equations 




then we get 



d G 

— (. x ) — lim — 
dx 5 a * 



y / rx+bx 

im — ( / 

^0 8x \J a 



fix) dx 



f 



fix ) dx 



= lim — 

8a— >0 OX 



^ rx+bx 



l 



fix) dx. 



The expression 



/»x+8x 

J x 



fix) dx 



(5.6) 



represents the area in the little strip between x and x + 8x (see Figure 5.2). Since 
fix) ~ fix) for this range of x the value of (5.6) is roughly Sx fix), and so we 
get 

G\x) ^ lim J- 8x fix) = fix): 

8a— >0 8x 



in other words G(x) is an anti-derivative of fix). This argument can be made 
precise if / is continuous (see Exercise 5.9). 

We now show how we can use an anti-derivative in order to calculate a definite 
integral (between two fixed limits) as in (5.3). If F is any anti-derivative of / then 
(d/dx)(E — G) = F' — G' = 0 and so F and G can only differ by a constant. 



Fix) = G(x) + c. 





5.2 General solutions and initial conditions 



25 



Since G(a ) = 0 (from its definition in (5.2)) we have F(a) = c, and so 

b 

f(x)dx = G(b ) 

= F(b) - c 
= F{b ) - F{a), 

which is (5.3). □ 

Because of the relationship between anti-derivatives and integrals, the notation 

I fix) dx, (5.7) 

(note the lack of limits on the integral) is often used as a shorthand to mean ‘an 
anti-derivative of /’. We will use this notation at times, but when we need to be 
more careful we will explicitly use a particular choice of anti-derivative F (x ) . 



5.2 General solutions and initial conditions 

Now let us return to our simple differential equation 2 

d y 

— = f(x). (5.8) 

dx 

Any anti-derivative F of / is a solution of this equation (y(x) = Fix)): hence we 
could simply write 

y(x) = J fix) dx. 

If we choose one particular anti-derivative F, then we know that not only is 
y(x) = F(x) a solution, but also y(x) = F(x) + c for any c. So as it stands (5.8) 
has many solutions. We say that 

y(x) — F{x) + c (5.9) 

is the general solution of the equation (5.8), since any possible solution of (5.8) 
can be obtained by choosing c appropriately. It should not be a surprise that there 
are many possible solutions; we can move a graph ‘up and down’ and not change 
its slope - all the curves in Figure 5.3, which differ by only a constant, have the 
same slope at any given x value. 

As illustrated in the figure, one way to pick out a particular solution is to give a 
point that must lie on the graph of y, in other words to specify the value y(jto) of 



2 



Usually we will not make explicit the dependence of dy /dx on x (and similarly for other derivatives). 
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5 ‘Trivial’ differential equations 




Fig. 5.3. Adjusting the constant c in (5.9) corresponds to moving the graph of F 
‘up and down’ and does not affect the slope. An initial condition (xo, >’o) will pick 
out one of the curves. 

y at some particular x value, xq. We refer to such a restriction 

y(* o) = yo 

as an initial condition , the idea being that we could construct the solution of (5.8) 
starting at (* 0 , >’ 0 ) and then drawing the graph by using the information about the 
derivative of y contained in (5.8). 

There are two ways to find the solution of (5.8) that satisfies y(x 0 ) = yo- You 
should make sure that you understand what follows, since we will use similar 
reasoning very often throughout the rest of the book. 

For the first method we do a little more than we have to: we find the general 
solution, and then solve a very simple algebraic equation to find the correct con- 
stant. We have seen that if we can find one anti-derivative F of / then the general 
solution of (5.8) is 



y(x) = F(x) + c. 

The particular solution that we want has v(xq) = yo, and so we need 
yo = y(-*o) = Fix 0 ) + c =>• c = yo - Fix 0 ). 

Thus the solution with y(x 0 ) = yo is 

yW = yo + Fix) - Fixo). (5.10) 

The alternative is to proceed more directly, and integrate both sides of 




5.2 General solutions and initial conditions 
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between xq and x. Then we get 



C x dv r x 

/ j-C’c) d.r = / f(x) dx, 



which gives 



yix) 



X px 

= / fix) dx. 

X=XQ J x 0 



Putting in the limits on the left-hand side this is 

>’(x) - y(x 0 ) = / fix) dx. 

JXQ 



(5.12) 



You should make sure that you are happy going straight from (5.1 1) to (5.12); we 
will generally skip the two intermediate steps. 

Since yix o) = yo, (5.12) gives the solution in the form 



yix) = yo+f fix) dx. (5.13) 



(The FTC shows that we do indeed have y'ix) = fix), and clearly yix o) = yo as 
required.) If we know that F is an anti-derivative of / then we can use the FTC, 



f 



fix) Ax = Fib) — Fia), 



(this was (5.3), and is just the usual rule for evaluating integrals) to write the solu- 
tion more explicitly as 



yix) = yo + Fix) - Fix o). 



Of course, this is the same expression that we obtained above in (5.10). Note, 
however, that in some cases the integral form in (5.13) may be the best that we can 
do, if it is not possible to find an explicit anti-derivative of /. 

We now look at some simple examples. 

Example 5.2 Find the general solution of the equation 

Ay 

— = x + lOsinv. (5.14) 

Ax 

What is the equation of the graph with slope x + 10 sin x passing through the point 

(77, 0)? 



In order to find the general solution we have to find an anti-derivative of 
v + 10 sin x. Using standard integrals this is ~ 10 cos x, and so the general 
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5 ‘Trivial’ differential equations 




Fig. 5.4. The graph of y(x) = \{x 2 — tt 2 ) — 10(1 + cos*). It passes through the 
initial condition (tt, 0), which is marked by a cross. 

solution of (5.14) is 

y(x) = \x 2 — lOcos* + c 

for any c. The one solution that passes through (tt, 0) must have 

0 — y(7T) = ^TT 2 10 + c =>• c=— ^TT 2 — 10, 

and so 

y(x) = \{x 2 — tt 2 ) — 10(1 + cos*). 

The graph of y against x is shown in Figure 5.4, along with the initial condition. 



Example 5.3 A curve passing through the point (1,0) has slope In*. What is the 
equation of the cun>e? 

We have to solve the equation 




Since we have already used the ‘long-winded’ method in the first example, let us 




5.3 Velocity, acceleration and Newton’s second law of motion 
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Fig. 5.5. The graph of y(x) — 1 + xlnx — x, with the initial condition (1,0) 
marked by a cross. 



do this one directly. We integrate both sides of the equation between 1 and x to 
give 



=i: 



y(x) — y(l) = / In xdx 



x In x — x 

Jx=l 

— (x ln.r — x) — (0 — 1) 



and so (since we want y(l) = 0) 

y(x) = l + xlnx — x. 

This solution is shown in Figure 5.5. □ 

In the next section we will give some more examples, this time more practically 
based, using Newton’s second law of motion (F = ma ). 



5.3 Velocity, acceleration and Newton’s second law of motion 

Newton formulated the calculus, and his theory of differential equations, in order 
to be able to write down and solve the mathematical models that resulted from 
his laws of motion. Since derivatives are essentially the ‘rate of change’, questions 
concerning velocities (the rate of change of position) and acceleration (the rate of 
change of velocity) are most naturally framed as differential equations. 
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5 ‘Trivial’ differential equations 



Newton’s second law of motion states that the change A p in the momentum p 
of an object is equal to F, the force applied, multiplied by the time At over which 
the force acts, 



A p = F At. 



Dividing by At and letting At tend to zero (this is, of course, a somewhat imprecise 
derivation) we obtain 



dp 
d t 



= F(t). 



Since the momentum is the mass m times the velocity v, i.e. p = mv, if the mass 
is constant we obtain 



dp 
d t 



d du 

= — (mv) — m — = F(t). 
dt dt 



The rate of change of v, dv/dt, is precisely what we mean by the acceleration, and 
so this equation is the familiar formula ‘F = ma’ written another way. 



Example 5.4 A car of mass m is travelling at a speed vq when it suddenly has to 
brake; the brakes apply a constant force k until the car comes to rest. How long 
does it take the car to stop, and how far does it travel before it comes to rest? 



Using Newton’s second law we have 



since the force acts to oppose the motion of the car. Rewriting this as v = —k/m 
and integrating both sides between times 0 and t we get 



v(t) - Vo = 




or 



v(t) - vo = 



kt 
m ’ 



and so 



v(t) = v 0 - 



kt 

m 



see Figure 5.6. The car stops when v(t c ) = 0; this implies that t c = mvo/k. 

Since the velocity v is the time derivative of the position x, v = x and we have 



dx 

dt 



kt 

= v(t) = v 0 . 



m 
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Fig. 5.6. The speed of a car that suddenly brakes: v (t) = vo — ( kt/m ), with 
m — 1000 kg and k = 6500 N and for initial speeds of 30 mph (solid line) and 
35 mph (dashed line). 




Fig. 5.7. The distanced moved by the car once the brakes have been applied. The 
choice of k and m is as for Figure 5.6, and again the solid line is for an initial 
speed of 30 mph and the dashed line is for an initial speed of 35 mph. 



Integrating both sides between t = 0 and t = t c (when the car stops) we get 



x(t c ) - x(0) = [ 

Jo 

Substituting for t c we have 



kt 

v 0 ) d t = 

m 



kt 2 

v 0 t - — 
2m 



f=0 



kt ^ 

= V 0 t c - — . 

2m 



x(t c ) - x(0) = 



mug 
2 ~k 



as shown in Figure 5.7. 
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5 ‘Trivial’ differential equations 



Since the stopping distance is proportional to the square of the speed, relatively 
small increases in speed will have a marked effect on the stopping distance. The 
stopping distance for a car travelling at 35 mph will be 49/36 of that for a car 
travelling at 30 mph, almost half as much again for just 5 mph extra speed, 3 see 
Figure 5.11. □ 



5.4 An equation that we cannot solve explicitly 

We remarked above that there arc many cases in which the best that we can do 
is to find the solution in the form of an integral, as in (5.13) where we wrote the 
solution of the general equation 

d y 

t - = /(*) w hh y(x 0 ) - }’o 

dx 



y(x) = yo+[ fix) dx. (5.15) 

Jx 0 

However, it can still be possible to describe qualitatively the behaviour of the 
solution. Here we consider a simple example, 

dx f 2 

— = e x(0)=xo. 

at 

Integrating both sides between times 0 and t gives the solution 



x(t) = xo + 




(5.16) 



This is as far as we can go without resorting to approximation, since there is no 
explicit form for the anti-derivative of e~‘ . 

However, it is known 4 that 



L 



oo 




d t = s/tt/2. 



3 Realistic values of m and k are m = 1000 kg and k = 6500 N (one newton is one kg m/s 2 ) which means that 
the stopping distances at 30 mph (« 13.4 m/s) and 35 mph (« 15.6 m/s) are 13.8 m and 18.7 m respectively (to 
fall prey completely to the British imperial/metric confusion, that is roughly 40 feet and 60 feet respectively). 
An extra 10 mph on the motorway means that the stopping distance at 80 mph is 64/49 of the stopping distance 
at 70 mph, over a quarter as much again: 225 feet (« 75 m) at 70 mph (^31 m/s) and 300 feet (^ 100 m) at 
80 mph (~ 36 m/s). 

4 This ‘Gaussian integral’ arises frequently, making it very frustrating that it cannot be evaluated explicitly. In 
particular the normal distribution, which is fundamental in the theory of statistics, is described by a bell-shaped 
curve whose equation is e - * / 2 /\/27r; statistical tables for the normal distribution are essentially based on 
evaluating the integral in (5.16) numerically. 
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Fig. 5.8. We cannot solve the equation explicitly, but we know that the solution 
always increases and tends to a value ^/tt/ 2 greater than its initial value. 

2 

So we can say, since x = e~'~ is always strictly greater than zero, that x(t) in- 
creases as t increases, and that 

x(t ) -h>- xq + Jk/2 

as t — > oo, see Figure 5.8. Even though we cannot write down an explicit form for 
the solution, we can still say exactly what happens ‘eventually’. In this way we 
can still understand something about the behaviour of the solution. This ‘eventual’ 
behaviour is often referred to as the long-time, or time asymptotic, behaviour. 



Exercises 



5.1 Find the general solution of the following differential equations, and in each case find 
the particular solution that passes through the origin. 

(i) 

d(9 

— = sinf + cost, 
dt 



(ii) 



dy _ 1 

dx x 2 — 1 



(use partial fractions) 
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5 ‘Trivial’ differential equations 



— At Inf, 



U 1 _ t 

— = e sin2f. 
df 

5.2 Find the function f(x) defined for —tt/2 < x < tt/2 whose graph passes through the 
point (0, 2) and has slope —tan x. 

5.3 Find the function g(x ) defined for x > — 1 that has slope ln( 1 + x) and passes through 
the origin. 

5.4 Find the solutions of the following equations satisfying the given initial conditions: 

(i) 

x — sec 2 t with x(tt/4) = 0, 



y' — x — jx 3 



with y(— 1) = 1, 



— =2sin“f with Qfn/A) = it/4, 
dt 



x — = l + x with V(l) = 1, 
dx 



[x(r)e-] 



with x(0) = 3, 



5.5 The Navier-Stokes equations that govern fluid flow were given as an example in Chap- 
ter 3 (see equations (3.1) and (3.2)). It is not possible to find explicit solutions of these 
equations in general. However, in certain cases the equations reduce to something 
much simpler. 

Suppose that a fluid is flowing down a pipe that has a circular cross-section of radius 
a. Assuming that the velocity V of the fluid depends only on its distance from the 
centre of the pipe, the equation satisfied by V is 



1 d / dV 



r dr V dr 



= -P, 



where P is a positive constant. 
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Multiply by r and integrate once to show that 

dV Pr c 

d 7 ~ 2 + r 

where c is an arbitrary constant. Integrate again to find an expression for the velocity, 
and then use the facts that (i) the velocity should be finite at all points in the pipe and 
(ii) that fluids ‘stick’ to boundaries (which means that V (a) — 0) to show that 

V(r)= j(a 2 -r 2 ), 

see Figure 5.9. (This is known as Poiseuille flow.) 

5.6 An apple of mass m falls from a height h above the ground. Neglecting air resistance 
its velocity satisfies 



m — = —mg v(0) = 0, 
df 

where v — y and y is the height above ground level. Show that the apple hits the 
ground when 




5.7 An artillery shell is fired from a gun, leaving the muzzle with velocity V. If the gun is at 
an angle 6 to the horizontal then the initial horizontal velocity is V cost), and the initial 
vertical velocity is V sin 6 (see Figure 5.10). The horizontal velocity remains constant, 
but the vertical velocity is affected by gravity, and obeys the equation u = —g. How far 




Fig. 5.10. Firing a shell at muzzle velocity V at an angle 6 to the horizontal. The 
shell follows a parabolic path. 





At 5mph over the 30mph limit, 

how much further does it take to stop? 




Twenty one feet, that's over six metres. 

At 30mph, you would have stopped in front of the first child. 

At 35mph, you would not only hit the first child, but you would still be 
doing 17mph when you hit the second child. 

It doesn't matter how good a driver you are. The faster you go, the longer 
it takes to stop. In this case the extra 5mph means the total stopping 
distance of an average saloon with ABS is an agonising ninety six feet, or 
twenty nine metres. Which brings you to the feet of the third child. 

30mph might feel slow, but only until someone steps out in front of you. 



Slow down. 



Fig. 5.11. A recent UK campaign to persuade drivers to cut their speed in town 
from 35 mph to 30 mph. The him at www.thinkroadsafety.gov.uk/ 
s lowdown/ download/ slowdown . mpg makes the point more forcefully. 
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does the shell travel before it hits the ground? (Give your answer in terms of V 
and 9.) 

5.8 In Dallas on 22 November 1963, President Kennedy was assassinated; by Lee 
Harvey Oswald if you do not believe any of the conspiracy theories. Oswald fired 
a Mannlicher-Carcano rifle from approximately 90 m away. The sight on Oswald’s 
rifle was less than ideal; if the bullet travelled in a straight line after leaving the rifle 
(at a velocity of roughly 700 m/s) then the sight aimed about 10 cm too high at a target 
90 m away. How much would the drop in the trajectory due to gravity compensate for 
this? (The initial vertical velocity v is zero, and satisfies the equation v — —g, while 
the horizontal velocity is constant if we neglect air resistance.) 

5.9 (T) This exercise fills in the gaps in the proof of the Fundamental Theorem of Calculus. 
Suppose that / is continuous at x, i.e. given any e > 0, there exists a S — 8(e) such 
that 



\x — x\ < S 



l/C*) - fix ) | < e. 



By writing 



j /-x+Sx 



= ±f 

Sx J x 



fix) = — / fix) dx 



show that for all S.v with \8x\ < <5(e) 



j rx+8x 



L 



fix) - Y~ / /(■*) 

ox 



< e. 



and hence that 



3J-S-0 &x 

You will need to use the fact that 

(*b 



j nx+8x 

lim — / fix)dx = fix). 



pb pb 

/ gix) dx < / I 
J a J a 



|g(x)| d.v < (b — a) max |g(x)|. 

xe[a,b] 
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Existence and uniqueness of solutions 



Because we arc going to spend some time trying to solve equations like 

dx 

— = f(x,t) (6.1) 

at 

we need to be sure that such equations will actually have solutions. Clearly it is a 
hopeless task to search for a solution of (6. 1) if the solution does not exist (hunting 
for a unicorn will take you a very long time). 



6.1 The case for an abstract result 

It is quite easy to write down equations that do not have any solutions, for example 

rt ^ d.X 

x~ + t ~ — = 0 *(0) = c 

d t 

does not have any solutions if c ^ 0: if t = 0 then the second term of the differen- 
tial equation disappears and we must have x (0) = 0. 

We have already seen that there arc many possible solutions when we want to 
find a function whose graph has a particular slope; the question of the uniqueness 
of solutions of a differential equation is somewhat subtle. However, we saw that by 
specifying a particular initial condition we could tie down one particular solution. 
So the problem that we will consider for our general theory will be the initial value 
problem (IVP), consisting of the differential equation supplemented by an initial 
condition, 

dx 

— = f(x,t) x (t 0 ) = x o. (6.2) 

dt 

If we suppose for a moment that the independent variable is time, then the re- 
quirement of uniqueness has a physical interpretation. Suppose that we specify an 
initial condition ‘now’, i.e. at time t = 0; then the existence of a unique solution 
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means that we can use the equation to predict the future, since the solution is 
uniquely determined for t > 0. In this context, uniqueness of solutions is equiva- 
lent to the requirement that our model be deterministic. 

Uniqueness is also useful since occasionally we may be able to guess what the 
solution of an equation is. If we substitute this guess in and it works, then it must in 
fact be the solution since we know that there is no other. We have already used this 
implicitly in Chapter 1 when we just checked that our solution N(t) = N s e k ^~ s) 
worked, and then assumed that it must be the only solution. 

As with existence, uniqueness is not automatic. The innocuous looking IVP 

Ax I At = sfx x(0) = 0 (6.3) 

has an infinite number of solutions. The ‘obvious’ solution is x{t) = 0 for all t > 0. 
But if you choose any value of c > 0, the function 

. [ 0 t < c 

xA,) = | « — c)*/4 t>c 

also satisfies the equation. Here the solution ‘waits around’ at ,r = 0, before even- 
tually ‘deciding’ (at time t = c) to wander off slowly to infinity. Some of the solu- 
tions of (6.3) arc shown in Figure 6.1. 

The issues of existence and uniqueness are real, and it is possible to come up 
with very simple equations in which they fail. The good news is that there is a very 
general theorem guaranteeing existence and uniqueness, with a hypothesis which 
is very simple to check. 




t 

Fig. 6.1. A number of solutions of equation (6.3). 
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6.2 The existence and uniqueness theorem 

The proof of the general existence and uniqueness theorem is beyond the scope 
of this book, and we will just state the result. However, if you arc interested 
Exercise 6.4 leads you through an outline version of the proof. 

In order to state the theorem properly we need to have a more precise idea of 
what we mean by a ‘solution’ of the initial value problem (6.2). The main point 
of the definition is that we allow for a solution to be defined only for some interval 
of t values, and do not require it to be defined for every t e M. 

Definition 6.1 Given an open interval I that contains to, a solution of the initial 
value problem 

dx 

— (t) = fix, t) with x(t 0 ) = xo (6.4) 

d t 

on I is a differentiable function x(t) defined on I, with x(to) = xo and x(t) = 
f{x, t)for all tel. 

The way that the definition specifies the interval on which the solution exists 
(rather than insisting that it be defined for every value of t e M) may seem pedantic 
at first, but we will soon see that this is necessary even for some very simple 
equations, since it is possible for the solution to ‘blow up’ in a finite time. 

But for now, given our formal definition of a solution, we can state the existence 
and uniqueness theorem. 1 

Theorem 6.2 If fix, t ) and df/dxlx, t ) are continuous for a < x < b and for 
c < t < d then for any xo £ la, b) and to e (c, d) the initial value problem (6.4) 
has a unique solution on some open interval I containing to- 

Essentially the result says that if the function f (x,t) is ‘sufficiently nice’ then 
the equation will have a unique solution, at least close to t = to (see Figure 6.2). 
However, the result tells us nothing about how large the interval is on which the 
solution can be defined. 

In almost all of the examples we meet, / will be ‘sufficiently nice’ ; but we have 
already seen one simple example in (6.3) for which there is no uniqueness. This 
does not contradict Theorem 6.2, since the derivative of x 1 / 2 is infinite at x = 0: 
when f(x) = x ^ 2 , we have fix) = fx” 1 / 2 , and this becomes infinite as x 0, 
so f is certainly not continuous at x = 0. 

1 In fact the conditions on / in the theorem are a little stronger than they need to be. It is only necessary that the 
function / is a Lipschitz continuous function of x, which means that 

\f(x,t) - f(y,t)\ < L\x -y\ (6.5) 

for some constant L. Any function with continuous first derivative is Lipschitz continuous (see Exercise 6.2), 
but not every function that is Lipschitz continuous has continuous first derivative (e.g. f(x ) = |x|). 




6.3 Maximal interval of existence 
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1 f * ) 

f o 

Fig. 6.2. Given an initial condition x (to) = xq, the existence and uniqueness the- 
orem only guarantees the existence of a solution defined on some open interval 
(marked by the bold line on the horizontal axis) containing the initial time to. 




Fig. 6.3. The derivative of x 1 / 2 , plotted here against x, is not continuous at zero, 
where it is infinite. 



6.3 Maximal interval of existence 

We now give an example showing that we need the freedom to specify the inter- 
val on which the solution of an equation exists if we want a result as general as 
Theorem 6.2: 



dx 
d t 



= x 



x(0) = xq. 



(6.6) 



Since x 2 and its derivative 2x are continuous, the equation certainly has a unique 
solution that exists in some open interval containing t — 0. We will see how to 
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6 Existence and uniqueness of solutions 



derive the solution of this equation in Chapter 8. For now observe that 

x(t) = (6.7) 

*0 - 1 

satisfies the equation (provided that xq ^ 0): clearly when t = 0 we have x(0) = 
vo, and differentiating gives 

r- |x -^ 1 -'f != (v ! q) = w ' )]2 ' 

so that the equation is satisfied. Since we know that the solution of the equation is 
unique, we must have the solution. 

Our solution has some interesting properties. If xo > 0 then the denominator 
is initially positive (at t = 0), but decreases as t increases until it reaches zero at 
time t — Xq 1 . This means that the solution, x(t). has become infinite by the time 
t = Xq X \ we say that it ‘blows up’ in a finite time. 

Things are much nicer, though, if we want to see where our solution came from 
in the past. We can decrease t (from zero) as much as we like, since as t de- 
creases the denominator becomes larger, and so the solution itself tends to zero 
as t —*■ — oo. So when xo > 0 we can define the solution of (6.6) on the interval 
(— oo, Xq - 1 ), but there is no way to define the solution on an interval that extends 
further into the future beyond the time t = x () 1 . We refer to (— oo, x^ 1 ) as the 
maximal interval of existence for (6.6). 

Note that Figure 6.4 also shows that solutions with xo < 0 tend to — oo as t 
decreases towards a finite t* < 0. When xo < 0 the maximal interval of existence 
is (x^ 1 , Too), and only for xo = 0 can we define a solution for all t e R (and then 
the solution is x{t) = 0). 

The two ill-behaved equations in this chapter (x = x 1//2 and x = x 2 ) should 
serve as cautionary examples as to the limitations of the existence and uniqueness 
theorem. That said, almost all the examples we meet in what follows will have 
unique solutions that exist at least for all t > 0. 



6.4 The Clay Mathematics Institute’s $1 000 000 question 

The questions of existence and uniqueness are somewhat dry, but there are still ex- 
tremely important mathematical models for which these issues arc not resolved. 
Outstandingly, it is still not known whether the Navier-Stokes equations that 




6.4 The Clay Mathematics Institute’s $1 000000 question 
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t 



Fig. 6.4. For positive initial conditions the solutions of i = x 2 blow up in a finite 
time, but exist for all negative values of r; while for negative initial conditions the 
solutions blow up for a finite value of t < 0 but exist for all t > 0. 



model the flow of fluids, 



P 



3u 

— + (u • V)u 
at 



— ix Au + Vp — f 



V • u = 0 



(cf. (3.3)), have unique solutions that exist for all positive times. 

These equations are the basis of computational design of everything that in- 
volves fluid flow; given that the term ‘fluid’ includes both liquids (in particular 
water) and gases (in particular air), numerical methods based on these equations 
are extremely important commercially. Clearly given the financial investment in- 
volved, people are confident that these equations really can predict the behaviour 
of physical systems, but currently we have no guarantee. Most tellingly, you can- 
not prove that a numerical approximation is ‘close’ to the ‘true solution’ if you do 
not even know that such a solution exists. 

For the year 2000 the Clay Mathematics Institute, based in America, announced 
seven Millennium Prize Problems; for the solution of any of these they will award 
a prize of one million dollars. 2 One of these problems is to determine whether 
or not the three-dimensional Navier-Stokes equations are indeed a good physical 
model, i.e. whether or not they have unique solutions valid for all positive times. 
There arc of course, two ways to win this prize: either to invent some insightful 
new mathematics that will prove the existence of unique solutions; or to dream up 
a single initial condition for which the solution breaks down. 



2 



See www . claymath . org/ index . htm 
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Exercises 

6. 1 Which of the following differential equations have unique solutions (at least on some 
small time interval) for any non-negative initial condition (x(0) >0)? 

(i) x = x(\ — x 2 ) 

(ii) x — x 3 

(iii) x = x 1 / 3 

(iv) x — x 1/,2 (l + x) 2 

(v) i = (1 + .r) 3 / 2 . 

6.2 (T) The Mean Value Theorem says that if / is differentiable on an interval [a, b] then 
f(a) — f(b) — (b — a)f'{c) for some c e (a, b). Suppose that f(x) is differentiable 
with \f'(x)\ < L for a < x < b. Use the Mean Value Theorem to show that for a < 
x, y < b we have 

\f(x)~ f(y)\<L\x-y\. 

6.3 (T) This exercise gives a simple proof of the uniqueness of solutions of 

x = f(x,t) x(to) = xo, (E6.1) 

under the assumption that 

\f(x,t) — f(y,t)\ < L\x — y\. (E6.2) 

Suppose that x(t) and y(t) are two solutions of (E6.1). Write down the differential 
equation satisfied by zit) — x(t) — y(t ), and hence show that 

j-|z| 2 = 2z[f(x(t), t) — f(y(t), f)]. 
at 

Now use (E6.2) to show that 

±\z\ 2 <2L\z\ 2 . 

at 

If dZ/dr < cZ it follows that Z(t) < Z(to)e c ^~ t0 ^ (see Exercise 9.7): use this to de- 
duce that the solution of (E6.1) is unique. Hint: any two solutions of (E6.1) agree when 
t = r 0 - 

6.4 (T) The proof of existence of solutions is much more involved than the proof of their 
uniqueness. We will consider here the slightly simpler case 

x — f(x) with x(0) — xo, (E6.3) 

assuming that 

I/W-/OOI <L\x-y\. (E6.4) 

The first step is to convert the differential equation into an integral equation that is 
easier to deal with: we integrate both sides of (E6.3) between times 0 and t to give 

x(t) = xo + [ f(x(t)) dr. 

Jo 



(E6.5) 
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This integral equation is equivalent to the original differential equation; any solution 
of (E6.5) will solve (E6.3), and vice versa. 

The idea behind the method is to use the right-hand side of (E6.5) as a means of 
refining any ‘guess’ of the solution x„(t) by replacing it with 



x n +i (0 = *0 + 

We start with xq ( t) = xq for all t, set 



[ f(x n 

Jo 






(E6.6) 



x\ (t) = xo + / /(xo)dr, 

Jo 

and continue in this way using (E6.6). The hope is that x„ ( t ) will converge to the 
solution of the differential equation as n —> oo. 

(i) Use (E6.4) to show that 



\x„+i (t) - x , 



n(t) | < L [ 

Jo 



I x„(t) - x n -i(t)\ df. 



and deduce that 



max \x n+ i(t) - x n (t)\ < - max \x n (t) - x„_i(f)|. (E6.7) 

rs[0,l/2L] 2 rs[0,l/2L] 



(ii) Using (E6.7) show that 



max |x, !+ i(f) - x n {t)\ < — — r max |*i(f) - xo(OI- 
re[0.l/2L] 2" _1 rs[0,l/2L] 



(iii) By writing 



x„(t ) = [. x n (t ) - x n -i (i)] + [x n -\ (t) - X n - 2 (?)] 
H 1- [x\ (0 - xo(i)] + xq it) 



deduce that 

1 

max \x n (t) - x m (t)\ < max \x\(t) - x 0 (i)| (E6.8) 

re[0, 1/2L] 2 N ~ 2 fs[0,l/2L] 



for all n,m > N. 

It follows that x n (t) converges to some function a'oo (t) as n — >• oo, and therefore taking 
limits in both sides of (E6.6) implies that 



*oo(7) = *0 + 




(t))dt. 



Thus Xoc (t) satisfies (E6.5), and so is a solution of the differential equation. The pre- 
vious exercise shows that this solution is unique. 
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Scalar autonomous ODEs 



For the most part when considering first order equations we will concentrate on 
finding explicit solutions. However, in this chapter we will see how, for the partic- 
ular class of equations of the form 



dx 
d t 



= fix). 



we can understand the solutions ‘qualitatively’, even if we cannot (or do not) write 
down their solutions explicitly. 

What this means is that instead of writing ‘x(f) = something’ we describe how 
the solutions behave, e.g. ‘any solution starting with x(0) between zero and one 
tends to v = 1 as t — > oo’ or ‘the point x — — \ is stable’. There is a very simple 
way to represent all this information about solutions pictorially, and the method 
essentially reduces to sketching the graph of the function / (in fact we only need 
to know where / is positive and negative). Nevertheless, the qualitative results we 
will obtain are completely rigorous. 



7.1 The qualitative approach 

The key observation is that the existence and uniqueness result of Theorem 6.2 
tells us that, provided / is ‘nice’, a solution of 

dx 

— = fix) (7.1) 

dt 

is completely determined by its value at any time t . The equation itself can then be 
used to determine whether x(t) is increasing or decreasing, depending on the sign 
of/. 

The easiest way to think of this kind of equation is to imagine that x(t) repre- 
sents the position of a particle moving on a line at time t. We can then talk about the 
‘velocity of the particle’ rather than the more cumbersome ‘rate of change of x\ 
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Fig. 7.1. A sketch of the function / against x, and the phase diagram for the 
equation x = fix). On the phase diagram the stationary points are represented as 
crosses and the arrows indicate whether the solution is increasing or decreasing. 



Whatever the equation really represents, we can use the particle idea while solv- 
ing it, and then reinterpret our results for the original application when we have 
finished. 

In order to understand how solutions behave we first find all the values of x at 
which the particle does not move; this happens at the points x* where fix*) = 0. 
As is often the case with fundamental ideas, such points have many names; we will 
call them stationary points } 

In regions where fix) > 0 the solution x(t) is increasing, and so the particle 
is moving to the right; similarly wherever f(x) <0 the solution x(t) is decreas- 
ing, and the particle is moving to the left. Note that the particle cannot reverse 
the direction in which it is moving; if it were to do this then at some time t* 
it would have to be instantaneously at rest, so that x(t*) = fixft*)) = 0. But 
then xit*) would be a stationary point, and so the particle would not be able to 
move. 

The simplest way to present all this information is to draw a line representing the 
v coordinate. Stationary points are usually indicated by a cross (x); we then draw 
arrows on the line indicating the direction in which xit) is changing: if fix) > 0 
then the particle will move to the right and i f / (x ) < 0 then the particle will move 
to the left. If we sketch the graph of / then it is easy to see the regions in which / 
is positive and negative. An example is shown in Figure 7.1. 

The picture of the line, with the stationary points and the direction of travel of 
the solution indicated, is known as the ‘phase diagram’ or ‘phase portrait’, which 
is shown on its own in Figure 7.2. (Figure 7.1 has two components; the phase 
diagram, and the sketch of / that makes it easier to draw.) With this kind of picture 



1 Other common terms are equilibrium points, fixed points and critical points. 'Equilibrium point' has a more 
physical flavour than the general tone of this book, and we will reserve the term ‘fixed point’ for use with 
iterated maps in Chapter 23 (a fixed point will be a point for which x* = f(x*)). We will use the term ‘critical 
point’ for a point at which a function F(x) has all its partial derivatives zero, see Appendix C. 
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{> X « X <3 X > X — <3 

Fig. 7.2. The phase diagram from Figure 7.1. 



qualitative behaviour of the solutions at a glance, even when we cannot write down 
the solutions explicitly. 



7.2 Stability, instability and bifurcation 

Looking at the phase diagram for the above example, we can see that some sta- 
tionary points are ‘attracting’ (nearby solutions approach), while some appeal - to 
be ‘repelling’ (nearby solutions move away). These ideas can be made mathemat- 
ically precise and are extremely important in applications. 

A stationary point is stable if when you start close enough to it you stay close 
to it. More precisely, a stationary point x* is stable if given any e > 0 there exists 
a 8 > 0 such that 



|je 0 -x*| <8 => \x(t)-x*\<e for all f>0. (7.2) 

v ' v " 

start close enough stay close 

The stationary points in the example above have the stronger property of being 
attracting , which means that if you start close enough you actually tend to the 
stationary point: there exists a S > 0 such that 



|*o _ x*| < <$ => x(t) — ► x* as t — ► +oo . (7.3) 

" v ^ v ' 

start close enough tend to 

In one-dimensional systems attracting points are stable (see Exercise 7.9) but 
this is not true in general. Conversely, an example with many stable stationary 
points that are not attracting is 



dx 
d t 



—x x < 0 

0 0 < x < 1 

1 — x x > 1 . 



Here, all the points in the interval [0, 1] are stationary points, and they are all 
stable. However, none of them are attracting, since there are nearby points that 
move no closer. The phase diagram is shown in Figure 7.3. 

A stationary point is unstable if it is not stable; this means that no matter how 
close (5) you start you will always move some stationary distance (e) away: there 
exists an e > 0 such that whatever S > 0 you take there is a point with |xo — x*| 
< 8 but | x(t) — x*| > e for some t > 0. 

The stationary points from Figure 7.2 are labelled according to their stability or 
instability in Figure 7.4. 




7.3 Analytic conditions for stability and instability 
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Fig. 7.3. The thick line consists entirely of stationary points, all of which are 
stable but none of which are attracting. 



E> — X < X — < X E> — X — <3 

s u u s 

Fig. 7.4. The phase diagram from Figure 7.2, but with the stationary points la- 
belled according to their stability type, S for stable and U for unstable. 





Fig. 7.5. The stability or instability of a stationary point x* can be determined 
from the value of fix*) provided that f'(x*) ^ 0. If f'(x*) < 0 then the sta- 
tionary point is stable, and if fix*) > 0 the stationary point is unstable. 



7.3 Analytic conditions for stability and instability 

There are very simple conditions on the derivative of / which will let us know 
whether a stationary point x* is stable or unstable without having to sketch the 
graph of /. 

If the graph of / near x* looks as shown in the left-hand side of Figure 7.5, 
i.e. if f'(x*) < 0, then the point will be stable, while if the graph of / looks 
like the right-hand side, i.e. if f'{x*) > 0, then the point will be unstable. Only 
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unstable/semi-stable 



< ^ X ^ < 



unstable/semi-stable 





Fig. 7.6. The various stability possibilities when f'(x*) = 0. 



when f(x*) = 0 is there any ambiguity; there are four possibilities, pictured in 
Figure 7.6. The top right and bottom left cases are sometimes called ‘semi-stable’, 
since the stationary point is stable ‘on one side’ and unstable on the other. 



7.4 Structural stability and bifurcations 

Observe that if fix*) ^ 0 then making small changes to the function / will not 
have a significant effect on the graph of / near x*. (The easiest ‘small change’ to 
imagine is adding or subtracting a constant, which will pull the graph of / up or 
down.) After the change there will still be a stationary point close to x* with the 
same stability properties (see Exercise 7.1 1). 

However, if fix*) — 0 then we can make small changes to fix) and drasti- 
cally affect the ‘picture’ near x*. For example, in the top right case of Figure 7.6, 
increasing fix) by any constant c > 0 will mean we no longer have a stationary 
point. 

When we make a small change to / but the phase diagram changes drastically 
we say that the equation has undergone a bifurcation. In these simple examples, we 
cannot have a bifurcation near x* unless fix*) = 0. When small changes to fix) 
cannot effect the qualitative nature of the phase diagram we say that x = fix) is 
structurally stable. 

We will look at a particular - example of a bifurcation in Section 7.6; watch for a 
stationary point with fix*) = 0. 



7.5 Some examples 

We now consider various examples using this graphical method. 




7.5 Some examples 
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Fig. 7.7. The phase diagram for the population model (7.4). Note that this picture 
includes solutions with p < 0, which, while mathematically sensible, are irrele- 
vant for this application. 



The equation 



7.5.1 A population model 



— = kp ( 1 — — ) with k.M > 0 

At V Ml 



(7.4) 



is a model for the change in the size of a population. 2 We will study this model 
in more detail in the next chapter, but for now we try to understand its qualitative 
behaviour. 

The first step is to find the stationary points. These occur where the right-hand 
side is zero, i.e. when 




= 0 , 



so they are p = 0 and p = M. If we sketch the graph of f(p) = kp( 1 — ( p/M )) 
then it is easy to draw the phase diagram, remembering that solutions move to the 
right whenever / > 0 and to the left whenever / < 0. The phase diagram is shown 
in Figure 7.7. It is easy to see from the diagram that provided that we start with a 
positive population then it will eventually settle down to the value at the stationary 
point p = M\ smaller populations will tend to increase, while hu ger populations 
will shrink towards this value. 

We can check to stability of the stationary points analytically by looking at the 
derivative of /, 



f\p ) = k - (2 kp/M). 



At the origin /'( 0) = k > 0, and the origin is unstable (as we expected), and 
f'(M) = — k < 0, confirming that the stationary point at p = M is stable. 

Note that this kind of solution can tell us what happens eventually, and how 
the population changes qualitatively. But it only allows us to make a very limited 



- There are some very reasonable objections to this model; in particular, unlike the number of people in a popu- 
lation, the variable p does not have to be an integer. We can get round this to some extent by claiming that p is 
‘the population in millions’, and then p can be a decimal. There are still values that do not correspond to whole 
numbers of people, but the equation will now be a good approximation. 
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t 



Fig. 7.8. The phase diagram (on the right) reflects the behaviour of solutions (a 
collection of which are shown on the left). 



number of quantitative predictions. For example, if we know that equation (7.4) is 
the right model for our population, but do not know the values of k and M. there 
is no way that we can use our phase diagram to find k and M given a collection 
of data. However, this is possible using the explicit solution, as we will see in the 
next chapter. 

Figure 7.8, shows how the phase diagram (rotated to be vertical in the figure) 
reflects the behaviour of the solutions, a collection of which arc plotted against t. 



7.5.2 Terminal velocity 

Sometimes we do not need an explicit solution to find the quantitative information 
we require. Here we will use the phase diagram to find the terminal velocity of a 
falling object. 

Suppose that a body of mass m is falling under gravity g and is subject to an 
air resistance proportional to the square of its velocity, kvr. The equation for the 
downward velocity v is 



m — = mg — kv\v\, (7.5) 

d t 

since gravity serves to accelerate the particle, and the air resistance acts in the 
opposite direction to v. 

Provided that the particle is moving downwards, so that v > 0, equation (7.5) 
becomes 




7.5 Some examples 
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Fig. 7.9. You can simply read off the terminal velocity v = y/mg/k. 



We can rewrite this as 



dr> 
d t 



k , 

= f(v) = g v 2 . 

m 



For this equation there is only one stationary point, when 




= 0 , 



i.e. when v = v* = y/mg/k. This point is stable, as can be seen by looking at 
the derivative of / at v *: since f'(v) = —2 kv/m we have f'(v*) = —lyfgkjm 
< 0 . 

The phase diagram is shown in Figure 7.9; it is clear that there is an attract- 
ing stationary point at v = Jmg/k. This is the terminal velocity, since v(t) ap- 
proaches this value whatever the initial condition. 

For a skydiver of mass 100 kg in freefall, we can take k ~ 1/3 kg/m and g ~ 
9.8 m/s 2 . It follows that the terminal velocity of the skydiver is 



v — VlOO x 9.8 x 3 ~ 54.2 m/s. 



7.5.3 What have we lost? 

We will now see that we are missing some, at times vital, information if we only 
rely on the phase diagram. The phase diagrams for the two equations x — \x\ and 
x = x 2 arc shown in Figure 7.10: although the equations arc different their phase 
diagrams arc the same. 

In the next chapter we will see how to calculate the solutions of both of these 
equations. For now we will assume that we know what these solutions are; for an 
initial condition x(0) = xq > 0, the solution of x — |x| is 



x(t) = xoe r , 
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> X — > 

Fig. 7.10. Different equations can have the same phase diagram; on the left is the 
phase diagram for x = |x|, and on the right the phase diagram for x — x 2 . 

while the solution of x = x 2 is 

1 

x(t) = — , 

x 0 ~ f 

(cf. (6.7)). 

The first solution increases to infinity, but is defined for all t > 0 (in fact for 
all t e M). However, we have already used the equation x = x 2 to show that it is 
possible for solutions to blow up in finite time ( x{t ) — > oo as t — > x^ -1 ). This finite 
time blowup behaviour is not captured in any way by our phase diagram. 

So although the phase diagram gives the correct qualitative behaviour, we have 
lost all information on the rates at which things happen. When the equation exhibits 
blow up of solutions in a finite time, this is particularly unfortunate. 




7.6 The pitchfork bifurcation 

We now consider the equation 

x — x(k — x 2 ), (7.6) 

where k is a parameter. By varying k we can study a whole family of differen- 
tial equations. We will see that the qualitative behaviour of the solutions of (7.6) 
changes drastically as k passes through zero. 

When k < 0 there is only one stationary point, that at x = 0. If we write 
/(x ) = x(k — ,r 2 ) then fix) = k — 3x 2 , and at the origin we have /'( 0) = k. It 
follows that for k < 0 this stationary point is stable. The phase diagram is shown in 
Figure 7.11, together with the graph of /. 

When k = 0 there is still only the one stationary point at x = 0, although now 
/'( 0) = 0. In order to determine the stability of the origin when k = 0 we have to 
sketch the graph of / (x ) = — ,r 3 . It is then clear - that the origin is still stable, and 
that the phase diagram is the same for k = 0 as it was for k < 0. see Figure 7.12. 

Since /'( 0) = 0 when k — 0. there is the possibility of a bifurcation as k 
changes from negative to positive (see Section 7.4). When k > 0 there are two 




7.6 The pitchfork bifurcation 
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Fig. 7.11. The phase diagram when k < 0, and the graph of /. 




Fig. 7.12. The phase diagram when k — 0, and the graph of fix) — — x . 




Fig. 7.13. The phase diagram when k > 0, and the graph of /. 

new stationary points at x = fz^/k. While the origin is no longer stable, since 
/'( 0) = k > 0, the new fixed points are both stable, since f'(±*/k) = —2k < 0. 
The phase diagram for this case is that shown in Figure 7. 13. 

You can see that the phase diagram has changed drastically as k has gone from 
being negative to positive. We have gone from having one stable stationary point 
for k < 0 to having three stationary points when k > 0; two of these arc stable, 
and the origin has become unstable. 

We can draw a ‘bifurcation diagram’ to show these changes. The idea is to draw 
a graph where the horizontal axis represents the parameter k, and for each value 
of k we plot the location of the stationary points on the vertical axis, using a solid 
line when they are stable and a dashed line when they arc unstable. This gives 
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Fig. 7.14. The pitchfork bifurcation; stationary points are plotted against k. Solid 
lines indicate stable points, and the dashed line an unstable point. 

the picture in Figure 7.14. For fairly obvious reasons this is known as a ‘pitchfork 
bifurcation’. 



7.7 Dynamical systems 

The qualitative approach we have adopted here is the main viewpoint used in the 
general theory of dynamical systems. A dynamical system has two components: the 
phase space (or ‘state space’), which consists of all possible ‘states’ of the system 
(for the scalar equations of this chapter this is the line M, covering all possible 
values of x), and the ‘dynamics’ which describe how these states change in time 
(for us the dynamics were determined by the solutions of the differential equation 
x = fix)). 

With the advent of more powerful computers there have been major advances in 
the theory of dynamical systems in recent years, and the subject received a lot of 
attention in the 1980s under the media-friendly ‘chaos’ banner. We will see more 
examples of dynamical systems later in the book. 

Exercises 

7.1 For each of the following differential equations draw the phase diagram, labelling the 

stationary points as stable or unstable. 

(i) x = — x + 1 

(ii) x — x (2 — x) 

(iii) x — (1 + x)(2 — x) sin x 

(iv) x — —x(l — x)(2 — x) 

(v) x — x 2 — x 4 





Exercises 
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7.2 For the equations in Exercise 7.1 determine the stability of the stationary points ana- 
lytically, by considering the sign of the derivative of the right-hand side. 

7.3 For all positive values of c find all the stationary points of 

d.x 

— = sin a- + c, 
d t 

and determine analytically which are stable and unstable. Draw the portion of the 
phase diagram between — tt and tt. There are three different cases, 0<c<l,c=l, 
and c > 1. You will need to be more careful with the case c — 1. 

7.4 A simple model of the spread of an infection in a population is 

H = -klH 
i = klH , 

where H{t) is the number of healthy people, / it) the number of infected people 
and k the rate of infection. Since (d/d t)(H + I ) = 0, it follows that the size of the 
population is constant, H + I — N, say. Substitute I — N — H in order to obtain a 
single equation for H(t), 

d H 

= -kH(N - H ). 

dr 

Determine the stability of the stationary points for this equation, and draw its phase 
diagram. Deduce that eventually all the population becomes infected. 

7.5 Consider the equation 

d.r t 

— = /(a) = a ~ -k. 
at 

Draw the phase diagram for the three cases k < 0, k — 0 and k > 0, labelling the 
stationary points as stable or unstable in each case. Find the stability of the stationary 
points using an analytic method when k > 0. Show that /'( 0) = 0 when k = 0. Why 
is this significant? 

Draw the bifurcation diagram, with k on the horizontal axis and the fixed points 
plotted against k, indicating stable fixed points by a solid line and unstable fixed 
points by a dashed line. (This is known as a saddle node bifurcation.) 

7.6 Draw the phase diagram for the equation 

a = g(.r) = kx — x 2 

for k < 0, k = 0 and k > 0. Check the stability of the stationary points by considering 
g'{ a), and show that the two stationary points exchange stability as k passes through 
zero. Draw the bifurcation diagram for this transcritical bifurcation. 

7.7 One equation can exhibit a number of bifurcations. Find, depending on the values of 
k, all the stationary points of the equation 

a = h (a ) = —(1 + a)(x 2 — k) 
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and by considering h'ix) determine their stability. At which points, and for which 
values of k, are there possible bifurcations? 

Draw representative phase diagrams for the five distinct parameter ranges k < 0, 
k — 0, 0 < k < 1, k = 1 and k > 1, and then draw the bifurcation diagram. Identify 
the type of the two bifurcations. 

In the remaining exercises assume that / is a C 1 function, i.e. that both / and d f/dx are 
continuous functions. Note that such an / is smooth enough to guarantee that the equation 
x — f(x ) with x(/o) = xo has a unique solution. You may also assume that the solutions 
are defined for all t > 0. 

7.8 (T) Let x(t) be one solution of the differential equation 

x = f(x). 

Show that 

(i) if f (x(t*)) — 0 for some t* then x(t) — x(t*) for all t e M (the solution is con- 
stant, and x(t*) is a stationary point); and hence 

(ii) if fix {t*)) > 0 for some t* then f(x(t)) > 0 for all t e R (the solution can- 
not ‘reverse direction’). Hint: Use the Intermediate Value Theorem: if g is 
a continuous function with g(a) < 0 and gib) > 0 then there is a point c 
between a and b with g(c) = 0. 

Of course, a similar result to (ii) holds if f(x(t*)) < 0 for some t*. 

7.9 (T) Show that for autonomous scalar equations, if x* is attracting then it must also be 
stable. Hint: use (ii) above, 

7.10 (T) Suppose that x(t) is a solution of x = fix) that is moving to the right. Show that 
either x(t) — > +oo, or xit) — > x*, where x* is a stationary point. (Hint: If x(r) does 
not tend to infinity then it is increasing and bounded above, and so tends to a limit 
x*. Show that in this case we must have fix*) = 0.) A similar result holds if x(f) is 
moving to the left, with +oo replaced by — oo. 

7.11 (T) Suppose that x = fix) has a stable stationary point at xo, with fix o) < 0. Let 
g be another C 1 function. Use the following scalar version of the Implicit Function 
Theorem to show that for e sufficiently small the equation 

x = fix) + egix) 

has a unique stationary point near xo which is still stable. 

Theorem. Suppose that h (x , e), dh/dx, dh/de are all continuous functions of both 
x and e. Suppose also that hix o, 0) = 0 and dh/dxixo , 0) ^ 0. Then there is an open 
interval I that contains xo such that for each e sufficiently small there is a unique 
solution v(e) e / of 

Hy(e),e) = 0, 



and yie) depends continuously on e. 




8 

Separable equations 



We now begin our survey of the various different classes of equations that we can 
solve explicitly. Both the ‘trivial’ equations 

dx 

— (0 = /(0 
at 

of Chapter 5 and the autonomous equations 

d.x 

- 7 ~(t) = f(x) 
at 

of the previous chapter arc particular cases of the separable equation 

7 ^ = fix) g{t) ( 8 . 1 ) 

at 

which we study in this chapter. 



8.1 The solution ‘recipe’ 

If you have already seen these equations, then you will probably be used to solving 
them in the following way. If these equations are new to you, take careful note; this 
is the practical way of finding a solution. However, there arc steps here that should 
make you uneasy. 

We start with the equation 

dx 

-77 = fix) gft). 
at 

Now divide by fix) and ‘multiply up by d t' to obtain 

1 

dx = git)dt. 

fix) 
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This is ‘separating the variables’, since we now have all the vs on one side and all 
the ts on the other. For the general solution we integrate both sides to get 

/7F) 4l = / s< ' )d '- (8 - 2) 

Alternatively, if we want to take into account an initial condition x (to) = xq then 
we integrate between the limits that correspond to times to and t : for the left-hand 
side these are x(to) and x(t), while on the right-hand side they arc just to and t. 
This gives 

rx(t) i ft 

/ —dx= g(t)di. (8.3) 

Jx o fix) J, 0 

We now use this recipe to find the solution of the equation we used in Chapter 6 
to show that the solutions of a differential equation can blow up in a finite amount 
of time (i.e. jt(f) — > +oo as t — »• t* < oo). At the time we had no method for 
solving this equation, and just wrote down the solution, but now we can use the 
separation method to find it for ourselves. 

Example 8.1 Find the solution of the initial value problem 

dx 0 

— = x x(0) = xo- 

dt 



If xo = 0 then x(t) = 0 for all t. Otherwise we can separate the variables to give 

1 

— dx = dt. 

x L 

Integrating between limits corresponding to times 0 and t, 

*x(t) i 



rx(t) ^ 

/ — r dx = / dt, 

Jx o * Jo 



we obtain 



Therefore 



1 

x 



1 x(t) 



= t. 






1 1 

H — t 



x(t) xo 



which simplifies to give 



x(f) = — 



X 0 -t 



as we claimed before (cf. (6.7)). 



□ 




8.2 The linear equation x — Ax 
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8.2 The linear equation x = Ax 

We now find the solution of the simplest possible linear differential equation. 




(8.4) 



with the initial condition x(to) = xq. 

This example is absolutely fundamental (the reasons for this will become appar- 
ent later) and you should really only have to solve this equation ‘long-hand’ once 
or twice before you arc happy to write down the solution with no calculation. 

First note that if xo = 0 then x(t ) — 0 for all t. Otherwise, if x ^ 0 then we can 
divide by x and ‘multiply up by dT to give 



dx 

— = A d t . 
x 



We now integrate both sides between the limits corresponding to the times to 
and t: that is, xo and x(t) on the left, and to and t on the right; and get 



which gives 1 



So we have 



r x(t) dx r ' _ 

/ V = 

Jx 0 Jto 



In lx I 



x(t) 



= Mt — to)- 






In |x(f)l — In |xo| = X(t — to), 
and taking exponentials (e to the power) of both sides gives 

Mill =c W-rn)_ 
kol 

To work out what to do about the modulus signs, the easiest thing is to draw the 
phase diagram. For the case A > 0 this is shown in Figure 8.1, from which we can 
see that x(t) and xo have the same sign. It follows that we can remove the modulus 
signs and multiply up to give 

x(t) =xoe k(t ~ ,o) . 



1 Many of the integrals in this chapter will involve logarithms, and the annoying modulus signs that come with 
them. We will have to take some care to work out how to remove them for our final answers; the most useful 
method is to use the phase diagram, as the examples show. 
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Fig. 8.1. The phase diagram for x — Xx when X > 0 (taking X < 0 would reverse 
the direction of the arrows). 



Note that the general solution of (8.4) is 

x(t) = Ae Xt , 



see also Exercise 8.5. 



8.2.1 Exponential decay and exponential growth 

We looked at the solution of equation (8.4) with X < 0 in Chapter 1, and applied 
it to the example of radioactive decay. We saw that the solutions decay to zero 
exponentially fast, and that the rate of decay could be characterised by the half- 
life; the solution halves in a fixed time. 

When X > 0 the solutions tend to infinity as t — oo, and increase exponentially 
fast. In this case the size of the solution will double after a fixed time, given by t 2 , 
where 



2 = e Xt2 , 

i.e. t 2 = In 2/a. In the following section we look at the use of this linear equation 
as a simple population model. 



8.3 Malthus’ population model 

The simple linear equation 

dp 

— = kp with k > 0 (8.5) 

d t 

was proposed in 1798 by the English economist Thomas Malthus as a basic model 
for population growth. Here the increase in the population is taken to be propor- 
tional to the total number of people, and k is a constant representing the rate of 
growth (the difference between the birthrate and the deathrate). This model pre- 
dicts exponential growth of the population, 

k(t-t 0 ) 



P(t) = p(to)et 




8.3 Malthus’ population model 
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so that its size grows without bound and will double every d years, where d = 
In 2 jk. This led Malthus to see war and famine as a possibly desirable check on 
this otherwise disastrous population explosion. 

We will now compare the predictions of this model with census data gathered 
over the last two hundred years. The population of Great Britain and Ireland in 
1801, 1851 and 1901 can be found in the results of the Census for each of those 
years: 

year population 
1801 16 345 646 

1851 27 533 755 ( ' ’ 

1901 41609 091 

We can use the data from 1801 and 1851 to estimate k. Our solution predicts 

/z (185 1) = p(mi)e 50k , 

and so 

10,(1851)-: In ,(1801) 

50 

This implies that the population will double roughly every 69 years (In 2 /k ~ 69). 

Using this value of k , our solution, illustrated in Figure 8.2, gives a reasonable 
prediction for the population in 1901: 

p ( 1 90 1 ) = p(1801)e 1(m « 46 million. 




Fig. 8.2. The UK population as predicted by Malthus’ linear model. The census 
values for 1801, 1851, 1901 and 2001 are indicated by crosses. 
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However, it vastly overestimates the population in 2001 as 

p(2001) = p(1801)e 2<m » 131 million, 

whereas the 2001 census found just below 59 million (in fact 2 the figure is the 
delightfully precise 58 789 194). To be consistent we should include the figures 
for the Republic of Ireland, since the data in (8.6) dates from before the partition 
of Ireland in 1921. The census held there in 2002 found a population of around 
4 million. 3 So the total figure for 2001 should be approximately 63 million. 
Malthus' model has predicted over twice this, so it turns out to be very unreli- 
able when we tty to extrapolate the population very far into the future. We will 
soon see another model that gives much more realistic results. 



8.4 Justifying the method 



We now give a careful justification of the ‘recipe’ we outlined in Section 8.1. In 
particular, you should have worried about the idea of ‘multiplying up by df since 
this kind of manipulation of infinitesimal quantities is extremely dubious. 

We start again with 

% = f(x)g«), (8-7) 

at 



and assume that fix) is sufficiently smooth to ensure that Theorem 6.2 guarantees 
the existence of a unique solution for any specified initial condition. 

First note that if x(t) is a solution of (8.7) with f(x(s )) = 0 for some s then in 
fact x(t) = x(s) for all t e M. This follows from the uniqueness of solutions; as- 
suming that x{t) = x(s) for all t implies that f(x(t )) = f (x (s ) ) = 0 for all (el, 
and so x(t) = 0 for all t, showing that this choice for x(t) solves the equation. 
Since solutions of the IVP are unique, this x(t) must be the only solution with the 
specified value of x (s ) . 

So either fix it)) — 0 for every value of t, or fix it)) ^ 0 for every value of t. 
We now treat the case f(x(t)) ^ 0 for all t, for which we can divide both sides of 
(8.7) by fix) to give 



1 dx 
f{x) d t 



= git)- 



( 8 . 8 ) 



Now, suppose that H (x ) is an anti-derivative of 1 /fix), i.e. 



H'ix) 



1 

fix)' 



2 See www. statistics . gov.uk/census2001/default . asp 

3 The exact figure was 3 917 336, see www. cso . ie / census /preliminary .details .html#pop . 




8.4 Justifying the method 
Then observe that by the chain rule (see Appendix C) 

d , d.r 1 dx 

and so (8.8) can be rewritten as 
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j~H(x(t)) = g(t). 
d t 

To find the solution we can integrate both sides with respect to t to give 

= J g(t)dt. (8.9) 

(To find x(t) explicitly we have to be able to invert H, i.e. solve the equation 
H(x) = z to obtain x in terms of z- In some cases the implicit form of (8.9) might 
be the best that we can do.) Since H is an anti-derivative of 1//, we could write 
this symbolically as 



S JiT) ix = S si,) 61 ' (810) 

which is precisely what we had before as equation (8.2). 

We now see how to recover (8.3) (see equation (8.11) below). If G is an anti- 
derivative of g then (8.10) reads 

H(x(t)) = G(t) + c, 

and so when we want to take into account an initial condition x (to) = *o we need 
H(x o) = G(to) + c =>• c = H(x o) - G(to), 
and the solution is 



H(x(t)) — H(xq) = G{t) — G(to). 



Using the method of evaluating an integral by anti-derivatives (which is formalised 
as (5.3) in the FTC) we can rewrite this as 




fix) 



dx 




git) dt, 



( 8 . 11 ) 



which agrees with the result of our more heuristic derivation above (equa- 
tion (8.3)). 
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8.5 A more realistic population model 

We now return to population modelling, but rather than allowing the unbounded 
exponential growth that resulted from Malthus’ model 




we impose a maximum sustainable size for the population. The idea is that any 
species (including ours) is limited by the availability of natural resources. We will 
find that this new model gives a much better estimate of the current population, 
even extrapolated from the century-old data we used above. 

The so-called ‘logistic equation’ is 



dp 
d t 




( 8 . 12 ) 



Interpreted as a population model, k is the growth rate of small populations; when 
p is small, p 2 is very small, so the equation is approximately dp/dt = kp, the 
model we had previously. The parameter M > 0 is the maximum sustainable pop- 
ulation; when p < M the population increases, and when p > M the population 
decreases. We drew the phase diagram for this equation in the previous chapter as 
Figure 7.7, and it will be useful to recall it now (see Figure 8.3) for use below. The 
phase diagram predicts that eventually the population will settle to its maximum 
sustainable level, p — M. 

We now solve the equation explicitly. Separating the variables gives 



M 

dp — d t, 

kp{M — p) 

where we have multiplied top and bottom of the left-hand side by M. Using the 
method of partial fractions on the left-hand side this becomes 



or 



1 

k 



1 1 ' 
p M - p 



dp = dt 



1 

- + 
P 




k dt. 



(8.13) 



X & X « 

0 M 

Fig. 8.3. The phase diagram for the population model (8.12). In line with the 
interpretation of p as the size of a population, only the values p > 0 are shown. 
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Since 

f 1 1 

/ - + — dp = In \p\ - In \M - p\ 

J p M - p 

We can integrate both sides of (8.13) between the limits corresponding to times to 
and t. 



L 



pit) i 
- + 



i 



dp = f k dt, 
P(t 0 ) P M ~P Jto 



to give 



- 


pit) 


- 


ln\p\ — In |M — p | 


= 


kt 


- 


p=pito) 


- 



jf=?o 



Putting in the limits of integration, 

In p(t) — In | M - p{t ) | - lnp(to) + In | AT - p(t 0 )\ = k(t - t 0 ) 



(since p(t) > 0 we do not need the modulus sign on In \p(t)\). Equivalently this is 



p(t)\M - pjt 0 )\ 

\M - p(t)\p(to) 



= k{t - to). 



From the phase diagram in Figure 8.3 it is clear that if pito) < M then pit) < 
M for all t, and similarly if pi to) > M then pit) > M for all t. So the sign of 
M — pit) does not change for each solution. It follows that we can remove the 
modulus signs, and then exponentiating both sides we obtain 



pjt)jM - pjt 0 )) _ 

Q^it—to) 

(M - pit)) pito) 



Finally, rearranging this gives 



Pit) = M 



pit 0 )e k(t ~ to) 

M - pito) + pito)z k{t ~ to) 



With some thought we can read from this explicit solution the same qualitative 
behaviour we see in the phase diagram. In particular, since c kl — oo as t — »• oo, 
we can deduce once again that pit) — > M as t — > oo. 

Since we have an explicit solution we can now estimate the parameters M and 
k that occur in the equation using the census data quoted above in (8.6). Once 
we know M and k we can then see what the quantitative predictions of the model 
are for 2001. The calculations to find M and k are just simple algebra, but arc 
not particularly instructive, so feel free to go straight to the values of M and k in 
equation (8.14). 
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Setting a = e 50k , po = p(1801), p\ — p(1851), and pi = p( 1901), our solu- 
tion requires 



P l 



MpQtx 

M - po + poa 



P2 = 



Mpoa 1 

M - po + poa 2 ' 



(Since there are now two parameters in our equation we need all the data from 
(8.6) to estimate them.) Rearranging both equations to find M in terms of a and 
equating we have 



popiia - 1) _ popiia 1 ~ 1) 

Poa - pi poa 2 - p 2 

and so 

Pl(pi-po) , P\i2P0Pi~ PiPi- PoP\) 

a — and M = -z 

poipi ~ P t) poP2 - Pi 



Using the correct values for po , p\ and pi gives 

a = e 50k = 2.0234 and M = 83.1 million 



(8.14) 



which implies that k ~ 0.014, similar to the value ( k ~ 0.010) we found for the 
simple exponential model (8.5). 

We can now use these values of k and M to predict the population in 2001. 
The value we obtain is 66.8 million, surprisingly close to the true figure (which, 
remember, is about 63 million); the solution is illustrated in Figure 8.4. Of course, 
there are good reasons for the discrepancy - among them two world wars and the 
invention of the contraceptive pill in the 1960s. 

Note that the constant M that arises in the model represents the maximum sus- 
tainable population; at 83 million this is still comfortably above its current level. 



8.6 Further examples 

We now treat some other examples. Note that often both the method of partial 
fractions and a quick sketch of the phase diagram arc useful tools. 



8.6.1 Partial fractions again 

We drew the phase diagram for the equation x = x(k — x 2 ) as part of our study of 
the pitchfork bifurcation in Section 7.6. Here we will consider the case k > 0, and 
so replace k by k 2 , which will make the algebra that is to come a little simpler. We 
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Fig. 8.4. Graph of the population of the UK and the Republic of Ireland, as pre- 
dicted by the logistic model using the census data for 1801, 1851 and 1901. The 
curve is our theoretical prediction, and the crosses show the exact values (we have 
chosen our parameters to ensure that the first three crosses lie on this curve). The 
dashed line is the maximum sustainable population predicted by our model. 



will solve the equation 



dx 
d t 



— x(k~ — x L ) 



with a general initial condition * (0) = xq. For the case k < 0 see Exercise 8.8. 
Separating the variables we have 



1 



dx — dt. 



x(k 2 — x 2 ) 

We can use the method of partial fractions to rewrite the left-hand side as 



1 



1 



1 



1 

- + 



1 



1 



x(k 2 — x 2 ) x(k — x)(k + x) K 2 \_x 2 (k — x) 2 (k + x) \ 

So we have 



1 



1 



1 

v ^ 2 (k — x) 2 (k + x) 



dx = k 1 dt. 



We can integrate this between the limits corresponding to times 0 and t to give 

1 1 



r x 

JXQ 



x(t) ^ 

- + 



2 (k — x) 2 (k + x) 



-L 



dx = k~ dt 
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+ — x « x ► x — « 

-k Ok 

Fig. 8.5. The phase diagram for x — x(k 2 — x 2 ). 



which is 



In \x | — i In \k — x \ — i In \k + x \ 



- x(t) 


to 

1 


JX=XQ 


L J 



t = 0 



or 



-\x(t) 



In 



\/\ K 2 ~ * 2 \ J.v=r 0 

Putting in the limits this becomes 



= K 2 t. 



In 



\x(t)\yj\K 2 — Xq\ 

\xqW\k 2 -x(t ) 2 1 
and exponentiating both sides we have 



= K 2 t, 



\x(t)\J\K 2 -X%\ _^ 2f 

\xo\V\k 2 ~ x(t ) 2 1 

Now if we square both sides and multiply up we have 

x(f) 2 k 2 - Xq | = xI\k 2 - x(t) 2 |e 2 ^ f . (8.15) 



We drew the phase diagram for this example in the previous chapter (see 
Figure 7. 13), and it is reproduced here as Figure 8.5. It is easy to see from the phase 
diagram that if Xo < K 2 then x(t) 2 < k 2 for all t and similarly if xo > then 
x(t) 2 > k 2 for all t. So we can remove the modulus signs and rearrange (8.15) to 
give 

x 2 ir 2 e 2l<2t 

x(t) 2 = = 

k 2 + XQ(e 2,t '“ r — 1) 

or 



x(t ) = ± 



1 + e 2l(lt {K 2 Xr ) 2 — 1) 



(8.16) 
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0 



t 

Fig. 8.6. Solutions of x — x(k 2 — x 2 ) on the left, and the corresponding phase 
diagram (rotated through 90 degrees) on the right. 

Whether we take the plus sign or the minus sign depends on the sign of the 
initial condition; for t — 0 we obtain 

*(0) = = ±|*ol 

and we have to choose the sign so that x(0) = xq. 

From the phase diagram we can see that for any xq > 0 the solution tends to k 
as t — oo, and for any xq < 0 the solution tends to —k; we can also recover this 
behaviour from our explicit solution, since if 0 < xo < k then e~' lK ~ t {k 2 Xq 2 — 1) 
is always positive, and decreases from its initial value to zero as t — > oo; it follows 
that x(t) increases from xq to k. Similarly if xo > k then e 2lc ‘ (k 2 x 0 2 — 1) is 
always negative, and increases up to zero as t — oo, so that x(t ) decreases to k as 
t oo. 

You can also see from the explicit solution that if |xo| > k then the solution 

r\ O 

will blow up as t \. t < 0, since then the expression (k~x 0 — 1) is negative (see 
Exercise 8.9 for more details). 

Figure 8.6 shows the solutions of the equation, along with the corresponding 
phase diagram, rotated to illustrate how the behaviour of the solutions matches the 
predictions of the phase diagram. 

8.6.2 Two competing species 

Later we will look at some simple models of competing species, and will come 
across such equations as 

dv y( 5x — 2) 

dx x (\ — 3y) 




72 



8 Separable equations 




Fig. 8.7. Curves on which In v + 2 In x — 3y — 5x is constant. 



This can be separated to give 



1 - 3>’ 5x — 2 
— ay = ax . 



Integrating both sides we have 



/(H d, -/( 5 - 



2 \ . 

- ) d.r; 
x 



taking x and y positive, because they represent the size of a population, we have 
no need of modulus signs in the logarithms arising from the integration. 

In y — 3y — 5x — 2 In x + c. 

We can do no better than this implicit solution relating x and y. However, we can 
represent the curves defined by 

F(x, y) = lny + 21nx — 3y — 5x — constant 

graphically, and these are shown in Figure 8.7. 



Exercises 

8.1 Solve the following equations: 

(i) x = r 3 (l — x) with x(0) = 3; 

(ii) y' = (1 + y 2 ) tanx with y (0) = 1; 

(iii) x = t 2 x (general solution); 

(iv) x = —x 2 (general solution); 
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2 9 

(v) for dy/d t = e ‘ y give the solution in terms of an integral and describe the 
behaviour of the solution as t -> +oo depending on the initial condition v (0) . 
You may assume that / 0 °° e -s ~ dj = 

8.2 Solve the linear equation 

x + px — C[ 

by separation of variables. 

8.3 Find the general solution of the equation 

xy — ky 



that is valid for x > 0. 

8.4 Find the function I (r) that satisfies 



d / 
dr 



= p(t)I. 



(Your answer will involve an integral.) 

8.5 Use the method of separation of variables to show that the general solution of the 
linear equation 



x = Xx 



is x(t) = Ae Xt for any /let 

8.6 In Exercise 5.7 we showed, neglecting air resistance, that an apple falling from a 
height h reaches the ground when t — *J2h/g. If we include air resistance then pro- 
vided that v < 0 the equation becomes 

dv 2 

m — = —mg + kv i/(0) = 0 
dr 

with k > 0. Show that 




and hence that the apple now takes a time 



r 



* 






to reach the ground. Check that this coincides with the answer with no air resistance 
(t* = y/ThJg) as A: — > 0. Hint: for small x, e v 1 + x and ln(l + x) x. 

8.7 Show that for k ^ 0 the solution of the differential equation 

d.r 2 

— = kx — x with x (0) = xo 
dr 



x(t) = 



k e k, xo 

xo(e kt — 1) + k 



is 




